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Abstract 

A  new  method  for  the  acceleration  of  linear  and  nonlinear  time  dependent  calculations  is 
presented.  It  is  based  on  the  Large  Discretization  Step  (LDS,  in  short)  approximation,  defined 
in  this  work,  which  employs  an  extended  system  of  low  accuracy  schemes  to  approximate  a 
high  accuracy  discrete  approximation  to  a  time  dependent  differential  operator.  Error  bounds 
on  such  approximations  are  derived.  These  approximations  are  efficiently  implemented  in 
the  LDS  methods  for  linear  and  nonlinear  hyperbolic  equations,  presented  here.  In  these 
algorithms  the  high  and  low  accuracy  schemes  are  interpreted  as  the  same  discretization  of  a 
time  dependent  operator  on  fine  and  coarse  grids,  respectively.  Thus,  a  system  of  correction 
terms  and  corresponding  equations  are  derived  and  solved  on  the  coarse  grid  to  yield  the  fine 
grid  accuracy.  These  terms  are  initialized  by  visiting  the  fine  grid  once  in  many  coarse  grid 
time  steps.  The  resulting  methods  are  very  general,  simple  to  implement  and  may  be  used  to 
accelerate  many  existing  time  marching  schemes. 

The  efficiency  of  the  LDS  algorithms  is  defined  cis  the  cost  of  the  computing  the  fine  grid 
solution  relative  to  the  cost  of  obtaining  the  same  accuracy  with  the  LDS  methods.  The  LDS 
methods  typical  efficiency  is  16  for  2D  problems  and  28  for  3D  problems  for  both  linear  and 
nonlinear  equations.  For  a  particularly  good  discretization  of  a  linear  equation  an  efficiency 
of  25  in  2D  and  66  in  3D  was  obtained. 
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1  Introduction 


In  recent  years,  interest  in  long-time  integration  of  partial  differential  equations  has  increased 
greatly  due  to  the  need  to  solve  diverse  problems  occurring  in  various  fields  of  science  and  engi¬ 
neering  such  as  fluid  mechanics,  aeroacoustics,  electromagnetism  and  others.  In  such  computations 
a  large  system  of  equations  has  to  be  evaluated  (explicit  schemes)  or  solved  (implicit  schemes) 
for  many  time  steps.  These  simulations  require  huge  computation  time  and  unless  more  efl&- 
cient  computational  methods  are  developed,  they  will  be  practically  intractable  in  the  foreseeable 
future. 

The  two  possible  approaches  to  employ  finite  difference  approximations  to  such  problems 
are  either  to  use  a  highly  accurate  scheme  (i.e.,  a  high  order  scheme  or  a  scheme  for  long-time 
integration  [10, 13])  on  a  grid  which  resolves  aU  the  physical  frequencies  occurring  in  the  problem, 
or  to  use  a  low  order  scheme  on  a  significantly  finer  grid,  or  a  combination  of  these  two.  The 
first  approach  seems  more  appealing  theoretically.  Indeed,  high  order  spatial  discretizations  [6, 
7,  13]  and  discretizations  for  long-time  integration  [10,  13],  as  well  ^s  high  order  time  marching 
schemes  [8]  have  been  in  the  focus  of  research  lately.  Although  significant  progress  has  been  made, 
the  two  main  problems  investigated  in  the  abovementioned  research  stiU  lack  general  solutions. 
The  appropriate  treatment  of  the  boundary  terms  in  high  order  Runge-Kutta  schemes  that  will 
maintain  the  interior  discretization  accuracy  still  requires  further  investigation  even  for  linear 
variable  coefficient  equations.  This  problem  is  mainly  of  a  theoretical  interest  as  it  has  only  a 
minor  effect  on  most  practical  computations  [8].  The  second  problem  is  the  lack  of  a  systematic 
method  for  constructing  numerical  boundary  conditions  of  a  required  accuracy  such  that  the 
resulting  discretization  is  time-stable.  This  is  a  major  obstacle  to  long-time  simulations.  The 
Large  Discretization  Step  (LDS)  methods,  presented  here,  offer  a  new  and  interesting  approach 
to  long-time  integration.  They  enable  to  obtain  a  fine  grid  accuracy  by  time  stepping  mainly 
on  a  coarse  grid  with  rare  visits  to  the  fine  grid,  at  a  cost  substantially  smaller  than  fine  grid 
simulation. 

In  some  cases,  the  huge  computational  cost  of  fine  grid  simulations  may  be  reduced  by  using 
such  a  grid  only  at  regions  where  it  is  required,  e.g.,  to  resolve  shocks,  and  employing  coarser  grids 
in  parts  of  the  computational  domain  where  the  solution  is  relatively  smooth  [1,  2, 3].  This  method 
of  local  mesh  refinement  for  systems  of  conservation  laws  has  been  reported  to  achieve  a  speedup 
of  up  to  a  factor  of  55  for  three  dimensional  problems,  relative  to  performing  the  computation 
on  a  uniform  grid  with  the  finest  mesh  employed  [1].  This  approach  assumes  the  scheme  has 
the  same  spatial  and  temporal  accuracy  and  does  not  seem  applicable  to  implicit  schemes.  The 
programming  effort  involved  in  generating  and  moving  the  fine  grid  patches  is  probably  the  cause 
for  the  limited  use  of  this  method. 

Multigrid  methods  have  been  employed  to  accelerate  time  dependent  computations  in  several 
ways.  The  naive  approach  is  to  use  an  efficient  multigrid  solver  for  implicit  time  marching  schemes. 
However,  in  this  approach  one  is  still  confined  to  the  fine  grid  time  step.  A  more  advanced  idea 
is  to  use  multigrid  in  time,  as  weU.  This  approach,  the  frozen  t  method,  was  successfully  applied 
to  parabolic  equations  [4,  5,  9].  There,  correction  terms  are  added  to  the  coarse  grid  equations 
such  that  one  can  time-step  on  the  coarse  grid  and  practically  obtain  the  fine  grid  solution.  This 
method  exploits  the  smoothness  of  the  change  in  the  solution,  typical  to  parabolic  equations, 
which  can  be  well  approximated  on  coarser  grids. 
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The  Large  Discretization  Step  (LDS)  methods,  first  introduced  by  the  authors  in  [11],  may  be 
viewed  as  a  generalization  of  the  frozen  r  method  aimed  at  accelerating  the  solution  of  hyperbolic 
as  well  as  parabolic  equations,  for  both  implicit  and  explicit  time  marching  schemes.  The  present 
work  investigates  the  LDS  approach  for  hyperbolic  equations;  a  class  of  equations  that  was  not 
previously  amenable  to  multigrid  methods.  Although  the  error  bounds  derived  in  this  work  apply 
to  any  time  dependent  equation,  the  algorithm  for  parabolic  equations  would  significantly  differ 
from  the  hyperbolic  solver.  Nevertheless,  it  is  expected  that  introducing  the  ideas  outlined  in  this 
work  to  parabolic  solvers  will  substantially  improve  their  performances,  as  well. 

The  present  paper  extends  the  preliminary  results  presented  in  [11]  in  several  important  as¬ 
pects.  The  algorithm  for  high  degree  LDS  was  significantly  improved;  resulting  in  a  more  efficient 
algorithm  that  requires  lower  order  intergrid  transfers.  The  method  was  extended  to  treat  non¬ 
linear  problems  with  the  same  efficiency.  Last,  problems  with  non-periodic  boundary  conditions 
were  solved. 

The  hyperbolic  LDS  identifies  two  grids,  a  coarse  representation  grid  on  which  aU  wavelengths 
occurring  in  the  physical  problem  can  be  well  resolved,  and  a  finer  computational  grid  which  is 
required  to  obtain  the  desired  accuracy  with  the  given  discretization  at  the  prescribed  final  time. 
The  LDS  method  employs  a  grid  possibly  finer  than  the  representation  grid,  yet  significantly 
coarser  than  the  computational  grid.  It  introduces  one  or  more  correction  terms  to  the  coarse 
grid  equations  and  a  system  of  equations  satisfied  by  these  terms  is  derived,  initialized  using 
the  fine  grid  and  solved  on  the  coarse  grid  to  yield  the  fine  grid  solution.  The  correction  terms 
are  integrated  on  the  coarse  grid,  hence,  their  accuracy  deteriorates  at  a  rate  determined  by  the 
coarse  grid  discretization.  However,  since  their  norm  is  significantly  smaller  than  the  solution 
norm,  they  can  be  effectively  used  for  many  coarse  grid  time  steps.  Thereafter,  the  fine  grid 
should  be  revisited  to  compute  new  initial  data  for  the  correction  equations. 

The  LDS  method  assumes  that  a  grid  which  resolves  all  the  physical  frequencies  occurring 
in  the  problem  as  well  as  a  discretization  suitable  for  a  fairly  long  simulation  time  are  given. 
However,  the  requirement  to  employ  the  same  discretization  for  substantially  longer  integration 
time  while  maintaining  a  desired  accuracy  necessitates  the  use  of  significantly  finer  grids.  The 
algorithm  solves  on  the  coarse  grid  an  extended  system  of  equations,  using  essentially  the  original 
time  marching  subroutines  (with  at  most  slight  modifications),  yielding  the  fine  grid  solution.  This 
programming  simplicity  renders  the  proposed  method  easily  applicable  to  any  problem  similar  to 
these  investigated  in  this  work,  provided  it  obeys  a  few  programming  conventions.  This  is  an 
important  feature  of  the  proposed  method. 

The  efficiency  of  the  LDS  is  defined  as  the  cost  of  computing  the  solution  on  the  fine  grid 
relative  to  the  cost  of  obtaining  the  same  solution  with  the  LDS  on  the  coarse  grid.  The  typical 
efficiency  achieved  in  this  work  was  16  for  2D  problems  and  28  for  3D  problems.  This  efficiency 
was  obtained  for  linear  problems  with  periodic  and  Dirichlet  boundary  conditions  and  for  the 
nonlinear  Euler  equation  in  a  periodic  domain.  A  particularly  good  discretization  yielded,  for  a 
linear  problem,  an  efficiency  of  25  in  2D  and  66  in  3D. 

The  organization  of  this  paper  is  as  follows.  Section  2  contains  a  heuristic  derivation  of  the 
method  for  both  linear  and  nonlinear  time  dependent  equations.  Section  3  presents  bounds  on 
the  error  in  the  LDS  approximation  for  linear  equations.  In  Section  4  it  is  shown  that  the  LDS 
approximation  maintains  the  stability  and  consistency  of  the  original  scheme.  The  LDS  algorithms 
are  described  in  Section  5.  In  Section  6  Fourier  analysis  is  employed  to  analyze  various  aspects  of 
the  algorithm,  in  particular,  to  obtain  the  necessary  orders  of  the  intergrid  transfers.  Section  6 
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presents  numerical  results  and  the  conclusions  are  given  in  Section  7. 


2  Heuristic  Derivation 


An  intuitive  and  informal  derivation  of  the  LDS  method  will  be  outlined  in  this  section.  A  rigorous 
derivation  for  linear  problems  is  given  in  Section  3. 

Consider  a  linear  time  dependent  system  of  equations  with  coefficients  possibly  dependent  on 
X  but  not  on  <, 


Tl  ^  ^ 

Mu  =  0 
u(a:,0)  =  uo(a^) 


for  X  e  f  €  [0,  T] 
for  X  G  do, 
for  X  € 


(2.1) 


where  SI  C  IR'',  and  ^  =  (,1^, sfj). 

Let  be  the  same  semi-discretization  of  equation  (2.1)  on  grids  h  and  JT,  respectively. 

Assume  that  the  fine  grid  is  required  to  obtained  the  desired  accuracy  at  time  T.  However,  instead 
of  solving  on  the  fine  grid 


(2.2) 


one  would  like  to  modify  the  coarse  grid  equation  such  that  it  will  yield  the  fine  grid  solution. 
The  coarse  grid  solution  satisfies, 

(2.3) 

A  correction  term  r  to  equation  (2.3)  is  sought  such  that, 

u'i  =  +  T  (2.4) 


where  denotes  a  restriction  of  the  fine  grid  solution  to  the  coarse  grid.  The  relative 

error,  ,  satisfies 

-  U^)t  =  L^(u^  -U^)  +  T  (2.5) 

Thus,  this  error  satisfies  the  same  equation  as  u^.  Moreover, 

(dt  -  L»)t  =  [dt  -  =  ri  (2.6) 


K  the  following  relation  holds,  which  is  reasonable  to  assume  when  well  approximates 


{dt  -  <  {dt  -  L")u^ 


(2.7) 


i.e.,  T\  <  r,  then  ri  may  be  neglected;  otherwise,  the  same  argument  implies  that  t\  satisfies  a 
similar  equation  to  resulting  in  a  larger  system  of  correcting  equations  (See  Section  3).  Thus, 
when  relation  (2.7)  holds  and  r  is  properly  initialized  the  system  of  equations 

u\  =  +  T  (2.8) 

.r  - 

Tt  =  L  T 
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yields  the  fine  grid  solution  on  the  coarse  grid  for  some  integration  time.  This  method  of  ex¬ 
panding  a  system  of  difference  equations  to  obtain  a  more  accurate  approximation  will  be  called 
an  LDS  approximation.  In  particular,  an  inflated  system  of  the  type  (2.8)  will  be  called  an  LDS 
approximation  in  Correction  Scheme  form. 

Introducing  the  new  variable  +  r,  the  system  (2.8)  can  be  written  as, 

-1-  (-2.9) 

uf  =  +  v^  - 

This  later  form  might  look  awkward;  however,  as  wiU  be  shown  next,  this  is  the  form  of  the  LDS 
for  nonlinear  problems.  This  representation  of  the  LDS  will  be  called  Full  Approximation  Scheme., 
following  the  multigrid  naming  conventions  [4]. 

In  Section  3  rigorous  error  bounds  on  such  approximations  for  linear  evolution  equations  are 
derived. 

For  nonlinear  problems,  only  a  heuristic  derivation  is  given.  Consider  the  nonlinear  evolution 
problem, 

Ut  =  P{x^u,—)  for  a:  €  te[0,T] 

M{u)  =  0  for  X  G  dft  (2.10) 

u{x^0)  =  uo{x)  for  a:  G 

where  D  C  and  ^  =  (gf^, 

Let  P^,P^  be  the  same  semi-discretization  of  equation  (2.10)  on  grids  h  and  H,  respectively. 
A  modification  of  the  coarse  grid  equation  is  sought  that  will  yield  on  that  grid  the  fine  grid 
solution,  i.e.,  a  forcing  term  is  required  which  wiU  satisfy 

u'l  =  P^{u^)  -t-  r  (2.11) 

where  denotes  a  restriction  of  the  fine  grid  solution  to  the  coarse  grid.  In  this  case 

the  relative  error  satisfies, 


=  P^{u^)-  P^{U^)  +  T  (2.12) 

w  P^{u^){u^  -U^)+T 

where  P^{u^)  is  a  linearization  of  P^  around  u^.  Assume  that  the  following  relation  holds,  which 
is  reasonable  if  P^  well  approximates  P^, 

{dt  -  P^{u^))  r  =  {dt-  Pf  (u'*))'  (u*  -  U^)  <  {dt  -  P^iu^))  {u^  -  U^)  (2.13) 

then  the  right  hand  side  of  the  r  equation  may  be  neglected. 

It  follows,  by  the  same  argument  as  in  the  linear  case,  that  the  system 

u'l  =  P^{u^)  +  T  (2.14) 

Tt  =  P^{u^)t 

yields  on  the  coarse  grid  the  fine  grid  solution,  for  some  integration  time. 
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This  formulation  of  the  IDS  is  inconvenient  to  use  as  it  requires  explicit  linearization.  Intro¬ 
ducing  the  new  variable  +  t  results  in  the  following  system  of  equations, 

=  P^{u^)  +  v^-u^  (2.15) 

v’}  =  P^{v^)  +  v^-u^ 

In  this  setting,  the  implementation  of  the  nonlinear  LDS  necessitates  only  minor  modifications  to 
the  original  time  marching  program. 

3  Approximation  Theorems 

The  LDS  approximation  introduced  in  Section  2  may  be  better  appreciated  once  the  initial  con¬ 
ditions  for  the  correction  equations  are  determined,  and  error  bounds  on  these  approximations 
are  derived.  These  two  issues  are  the  subject  of  the  present  section. 

First,  an  error  bound  is  obtained  for  a  semi-discrete  approximation  in  a  restricted  setting.  This 
restriction,  a  commutativity  assumption,  is  introduced  merely  to  maintain  a  simple  and  intuitive 
presentation.  Subsequently,  error  bounds  for  semi-discrete  and  fully  discrete  approximations  are 
derived  without  this  superfluous  assumption. 

The  analysis  will  be  performed  for  linear  equations  with  coefficients  which  may  depend  on  x 
but  not  on  f,  of  the  form, 

Q 

Uf  =  L{x,—)u+F{x)  for  a:  €  fi  t€[0,T] 

Mu  =  0  for  a:  €  dn  (3.1) 

u(a;,0)  =  uo{x)  for  a:  €  D 

where  D  C  and  ^  =  (gl^,  Jj, ...,  g|^). 

3.1  Semi-Discrete  Analysis 

3.1.1  Motivation 

Consider  a  stable  semi-discretization  of  a  linear  homogeneous  initial  boundary  value  problem  of 
the  form  (3.1)  given  by, 

^  =  0  (3.2) 

dt 

with  initial  conditions  u^(0)  ==  Uq. 

Let  be  an  approximation  to  e.g.,  a  coarse  grid  representation  of  the  fine  grid  operator. 
Define  the  system 


with  initial  data 

vl^{0)={L^-i^yu\0)  (3.4) 
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Henceforth,  an  approximation  of  a  system  of  equations  by  an  enlarged  system  of  the  form  (3.3) 
will  be  called  an  LDS  approximation  of  degree  m. 

Assume  that  and  commute.  The  solution  of  this  system  of  linear  ordinary  differential 
equations  is  where  denotes  the  above  system.  The  matrix  has 

a  block  Jordan  form;  hence,  an  explicit  expression  for  VQ^{t)  is  given  by. 


4,m(0  = 


{L^  - 


u^(0)  = 


(m+  1)! 


u*(0) 


for  some  ^  G  [0,t]. 

Assume  that  ||e^^*w^(0)||  <  C^e^^*l|u^(0)||,  for  constants  then 


-  <„(<)|| 


||(X*  -  ||i(.  _  £k||m+l,m+l 


(m  +  1)! 


{m  +  1)! 


&e^''%u^{<d)\\  (3.5) 


The  bound  (3.5)  implies  that  for  any  fixed  final  time  T,  lim„i_oo  ll'W^(r)  -  Vq^^{T)\\  =  0,  with 
convergence  rate  depending  on  the  magnitude  of  the  relative  truncation  error,  \\L^  -  L^\\. 

This  bound  also  suggests  that  when  the  relative  truncation  error  is  small,  then  there  exists  a 
time  To  such  that  the  error  in  the  approximation  of  a  fixed  degree  m  is  small  for  T  <Tq.  This 
observation  motivated  the  LDS  algorithm  and  enables  its  high  efficiency.  In  this  algorithm 
and  stand  for  the  same  discretization  on  the  fine  and  coarse  grids,  respectively.  The  algorithm 
computes  initial  conditions  for  the  correction  equations  using  the  fine  grid,  then  marches  with  the 
enlarged  system  on  the  coarse  grid  as  long  as  the  LDS  error  relative  to  the  fine  grid  solution  is  of 
the  same  magnitude  as  the  error  in  the  later  solution.  Then,  the  fine  grid  is  revisited  to  compute 
new  initial  data  for  these  equations.  In  this  manner  the  fine  grid  accuracy  can  be  obtained  when 
time  marching  mainly  on  the  coarse  grid. 

The  identity  in  (3.5)  can  be  used  to  obtain  the  following  inequality. 


„(()«< 


||e^'’‘||  IKZ,'-  -  X^)'“-»1|'-(0)||  ;”*+■ 

(m  +  1)! 


The  stability  of  the  semi-discretizations  considered  implies  that  there  exists  a  mesh  size  ho  and 
constants  C, /?  such  that  for  aU  grids  with  mesh  size  h  <  ho, 

||e^'‘^u"(0)||<Ce^'||u>)l|  (3.7) 

Therefore,  for  h  <hQ  the  inequality  (3.6)  implies  that, 

„„»(,) .  (3.8, 

yra  +  Ij!  ^  ^ 

Hence,  for  a  fixed  degree  m  and  fixed  integration  time  T,  the  error  in  the  LDS  approximation 
satisfies,  ||u^(f)  —  ^^0,771(011  =  0(h(”^'*‘^)^),  provided  {L^  —  L^)  —  0{hP).  Furthermore,  if  p  >  0, 
which  is  reasonable  to  assume,  this  bound  suggests  that  the  LDS  error  decreases  as  h(”^+^)p  as 
mesh  is  refined. 

In  the  next  section,  an  error  bound  is  proved  for  non  homogeneous  equations  without  the 
commutativity  assumption 
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3.1.2  Error  Bound 

Let  X*  be  a  discretization  of  the  spatial  operator  X  in  (3.1).  Let  be  an  approximation  to  L^, 
and  denote  -  A^.  Intuitively,  X^  may  be  viewed  as  a  fine  grid  discretization  and  A^ 

a  coarse  grid  approximation  to  X^.  However,  it  should  be  noted  that  all  the  operators 
are  defined  on  the  same  grid.  The  semi-discretization  of  equation  (3.1)  may  be  written  as 

uh  _  ^  ^  jh 

u^(0)  = 

Denote  the  solution  of  this  problem  by  U'''(t)  €  , where  N  is  the  number  of  grid  points. 

Assume  the  following  inequalities  hold. 


||(/'■(^)||  <  Ce'>'(||u‘||  +  ||/‘||) 


(3.10) 

(3.11) 


for  constants  C,/}  which  may  depend  on  h.  The  stability  of  the  semi-discretizations  considered 
implies  that  for  fine  enough  grids  (3.10)-(3.11)  may  be  bounded  independently  of  h. 

Define  the  system  of  equations 


l,7n 


(3.12) 


^oV(O) 

V;^-i,ni(0) 


(3.13) 


//  =  B^f  0<j<m 


with  the  Bj  defined  inductively  by 


(3.14) 


B^  =  I 

Bf^^  =  [Bf,A’^]ABfB>^  i>0 


(3.15) 


It  will  be  shown  that  the  first  component  of  the  solution  of  this  system,  approximates 

More  precisely,  a  bound  on  1|Io^to(0  ”  ^^^(011  presented  which,  for  fixed  t,  tends  to 


zero  as  m  00, 


The  vector  (jSq  C^^(t), . . . ,  B!^U^{t))  satisfies  the  equation 


B^U^ 


Bt-^U' 


A’^  I  0  . . . 
0  A^  I  ... 


B^U^ 


Bt-iV' 

BtU^ 


fm-\ 


(3.16) 
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(3.17) 


BLid^O) 


It  will  be  shown  that  converges  to  BjU’^itj  as  m  -r  od.  In  particular,  converges 

to  U^(t)  as  m  ^  00. 

Define  It  satisfies  the  following  equation 


_d 


■"m  — 1,77Z 
6^ 


I  0  . . . 

0  I  ... 


-"m— 1,771 
6^ 

'"771,771 


si+iU>^ 


(3.18) 


4-(o) 


d-l,m(0) 

em,m(0) 


The  solution  of  this  system  satisfies 


Xt)= 

J  0 


(3.19) 


(3.20) 


IM  =  f  0  <  j  <  m  -  1 

,/  0 


The  norm  ||€o^(t)||  is  the  sought  error  bound. 


|e‘ ,„(l)ll  <  ll<+,ll  j(V<'-)||(/‘(s)||<Is  <  ||Bi+,||  (||«5||  +  ll/'-ll)  C(c«' 


By  induction  one  obtains 


<  l|5‘+rll  (llaSlI  +  ll/'-ll)  C 


^771  — /+1 

(m  -  /  +  1)!* 


(3.21) 


(3.22) 


(3.23) 


The  following  theorem  was  proved. 

Theorem  1  Let  U^{t)  he  the  solution  of  (3.9)  and  (FqVCOi  •  •  -  i  K^,m(0)  that  of  (3.12)- 
(3.13)  Then 


lllSw)')  -  £''‘(011  <  lIBw+iit  (ll«Sll  +  ll/'•||)  C 


{m  +  1)! 


(3.24) 


A  different  bound  on  the  error  in  the  IDS  approximation  will  be  derived  next.  Denote 


^^(t)  =  sup 

0<5<t 

(3.25) 

*S.+i(<)  =  “P  lisi+iff'c*)!! 

0<5<t 

(3.26) 

The  exponent  /?  may  depend  on  the  mesh  size  h.  However,  for  a  grid  fine  enough  the  stability  of 
the  semi- discretization  implies  that  the  solution  can  be  bounded  independently  of  h.  Therefore, 
it  will  be  assumed  that  the  grid  is  fine  enough  for  this  property  to  hold.  Then, 

(3.27) 

and  . , , 

IMLWII  <  [«*(<)]  _;  +  l)! 

(3.28) 

Therefore, 

lkoV(i)ll  <  [«‘(<)] 

(3.29) 

This  bound  is  very  crude  and  can  be  easily  improved.  However,  it  implies  that  for  a  fixed  final 
time  T  and  degree  m,  if  Bm+i  =  0{hP^+^)  with  pm+i  >  0,  then  =  0{hP’”+^)  and  the 

error  in  the  LDS  approximation  decreases  like  as  mesh  is  refined. 

The  bound  (3.29),  suggests  that  looking  at  may  give  an  insight  into  the  accuracy  of  the 

approximation.  From  the  definition  of  Bj,  it  follows  that  B^  =  B^  and  B2  =  +  (5^)^. 

In  the  following  two  examples,  is  interpreted  as  a  coarse  grid  approximation  to  the  fine 
grid  spatial  discretization  A^  A  B^. 


Example  1  Consider  a  discretization  of  the  two  dimensional  wave  equation 

utt  =  Am 

when  transformed  into  the  system 


(3.30) 


(3.31) 


This  discretization  of  the  wave  equation  was  successfully  used  in  [12]  to  solve  problems  in  elasticity. 
Assume  that  the  Laplace  operator  is  approximated  by  a  second  order  scheme.  Then,  up  to  higher 
order  terms, 

(  A  +  a  +  dyyyy)  0  )  ’  ^  ^  +  Oyyyy)  0  ) 

where  a  denotes  a  generic  constant.  Since  =  0,  it  follows  that 


Bl  =  [5^A"]  = 


_  f  O  (h  H  )  {dxxxx  “k  ^yyyy  > 

0 


H^)  (dxxxx  +  dyyyy)  ^ 


It  can  be  seen  that  B2  consists  the  same  terms  as  the  relative  truncation  error.  Thus,  one  can  not 
expect  that  the  LDS  of  degree  one  will  yield  an  approximation  more  accurate  than  the  coarse  grid 
discretization.  Indeed,  applying  the  LDS  algorithm  to  this  discretization  shows  that  the  error  in 
the  LDS  of  degree  one  solution  grows  like  the  coarse  grid  error. 


Example  2  Consider  the  linearized  Euler  equation 


Pt  -  a-V  +  c{u^  +  Vy) 

ut  =  a  •  V  +  c  Pj  (3.34) 

Vf  =  a-V  +  cpy 


where  a  =  (01,02)  is  a  two  dimensional  vector.  Assume,  that  the  spatial  operator  is  discretized 
using  a  second  order  scheme.  Then,  up  to  higher  order  terms. 


Ah 


(  o  •  (V  +  a  {pxxx  +  ^yyy))  c{dx  + a  dxxx) 

c{d^  +  aH^d^^^)  a-{V  +  aH^d^^^  +  dyyy)) 

\  C{dy-\-<XH  dyyy)  0 


Cidy  +  aH^dyyy)  \ 

0 

o  •  (V  +  o ^  j 


=  {h?  -  H^)  ( 


^  dxxx  H"  ^  ^yyy) 
C  Ol  3xxx 

C  a  dyyy 


C  Oi  dxxx 

d  ‘  (^Cx  3 XXX  ^  ^yyy  ) 

0 


ca3. 


yyy 


d  *  (of  3xxx  ^  ^yyy^  J 


(3.35) 


where  a  is  a  generic  constant.  It  can  be  easily  seen  that  f  consists  of  sixth  order  mixed 
derivatives  and  fourth  order  powers  of  H,  h. 


/  0  0  0  \ 

=  0  0  p  \  (3.36) 

\  0  -p  0  J 

with 

Tj  —  C  a  [h  {3yyyx  ^  3xyyy  +  3yyyxxx  ”  3xXXyyy^  (3.3^) 

For  smooth  solutions, 

P=c‘^a^  {h^  -  H^fidyyy^  -  O^^yyy)  (3.38) 

Hence,  the  error  bound  (3.29)  implies  that  for  smooth  data  the  LDS  of  degree  one  yields 
a  significantly  smaller  error  than  the  coarse  grid  operator.  Indeed,  our  numerical  results  (see 
Section  7)  show  that  the  LDS  algorithm  for  this  equation  yields  the  fine  grid  solution  on  the 
coarse  grid. 


3.2  Fully  Discrete  Analysis 

Consider  a  stable  finite  difference  approximation  to  (3.1)  of  the  form 

t/"+i  =  (/  +  kA)U'^  +  kBU^  +  kf  ,  , 

(f»U„ 

where  k  —  /Si. 

Assume  that 


||C/”1|  <  e'^"^(||uotl+ ll/ll) 
\\I  +  kA\\  < 
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(3.40) 

(3.41) 


In  an  analogous  way  to  the  semi-discrete  case,  define  (Io%5  •  •  •  5  Kn.m)  the  solution  of 


yn-fl  \ 

(  I  +  kA  kl  0  . . . 

K"  \ 

(  kfo  ^ 

: 

0  I  +  kA  kl 

I 

1 

yn+l 
^  m~l,m 

m— l,m 

kfm-l 

y«+i  ) 

^  m,m  / 

y  I  +  kA  j 

i 

\  m,m  / 

\  k  fjYi  j 

(3.42) 


^  BqUo  ^ 

y^-l,ra 

— 

Vim 

\  -5m  '^O  / 

fj  =  Bjf 

0  <  j  <  rn 

with  Bj  defined  recursively  by 

Bo  =  I 

Bj+i  —  [BjTA]  +  BjB  j>0 

Consider  the  vector  (BoV^, Bm-iU'^,BmU'^)-  It  satisfies 

\  /  BqU^  \  ( 


(  \  (  I +  kA  kl  0 

0  I  +  kA  kl 


Bm-xU^^^ 

\  B^U-+^  ) 


lAkA  J 

with  the  same  initial  condition  as  {Vq  .^,  . . . , 


Bm-lU^ 

V  BmU-  ) 


(3.43) 


kfo 


(3.44) 


(3.45) 

(3.46) 


\ 


k  fm—l 
\  kfm  +  kBm+lU-)  ) 


(3.47) 


The  error  =  Bjlf^  -  satisfies 


/  eg+i  \  /  I  +  kA  kl  0 


\  eZ+l  } 

\  m,m  / 


0  /  +  kA  kl 


.  \ 


lAkA  j 


I 


^0,m 


pTl 

\  eZrr,  I 


0 

V  kBm+lU^  ) 


(3.48) 


^0,m  ^ 

0  N 

>0 

m— l,m 

0 

e“  J 

^m,m  / 

(3.49) 
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These  equations  give, 


=  Y.ii+kAr-^-^kBm^,u^ 

j=0 


(3.50) 


0  <  /  <  m  —  1 


(3.51) 


The  following  theorem  can  be  proved  by  bounding  the  solution  of  this  discrete  system. 
Theorem  2  Let  t/”  be  the  solution  of  (3.39)  and  {Vo^m{t), . .  .,Vm,m{t))  be  that  of  (3-42)- 
(3.43).  Then 

\\Vo%  -  U'W  <  (||f/»||  +  ll/ll)  (3.52) 

Proof  :  The  following,  more  general,  formula  will  be  proved 

Iktell  S  (l|£f°ll  +  ll/ll)  l|B««ll  ^  for  0  <  i  <  m  (3.53) 

Theorem  2  is  the  particular  case  of  /  =  0. 

The  proof  follows  by  induction  on  m  -  First  consider  the  case  m  -  /  =  0,  i.e.,  I  =  m. 
According  to  formula  (3.50) 


el^^  =  J^{I  +  kAr-^-^kB^+rU^ 

j=o 


(3.54) 


Therefore, 


\^m,m\\  < 


^  \\{i+kAr-^-^\k\\Bm+i\\  iit^i  <  (iif^°ii  +  ll/ll) 

j=0  j=0 

(||C^°II  +  ll/ll)  ||5„^+l||(n^r)e^("-l)‘  (3.5.5) 


Assume  the  bound  (3.53)  is  correct  for  m  ~  /  =  m  —  1;  that  is, 


<  (l|t^°ll  +  ll/ll)  \\Bm+l\\r]k^e^(--^^'^ 


(3.56) 


From  relation  (3.53)  it  follows  that. 


j=0 


Therefore,  by  the  induction  hypothesis  (3.56), 


K..II  <  (3.58) 


<  l|t^°ll  +  ll/ll  l|5: 


m  +  1 


J^m+l  g/?(n-rn-l)  k 


(3.57) 


(3.59) 
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The  last  inequality  follows  from  the  following  identity  which  can  be  proved  by  induction  on  n. 


n— 1 


su  = 


i=o 


n 

,m  +  1, 


The  relations  (3.50)-(3.51)  give  rise  to  a  somewhat  different  bound,  as  well.  Denote, 


(3.60) 


'$'*(n)  =  sup  11(7 -1- A: A)'^ II 

0<j<n 

(3.61) 

=  sup  ||B^+iC^^|| 

0<i<n 

(3.62) 

Then, 

lie”  <  ink)^\n)^Urin) 

(3.63) 

and 

\Km\\  <  {nkr^^  [^\n)]"  $^+i(n) 

(3.64) 

The  stability  of  the  discretization  implies  that  for  a  fixed  time  T,  if  the  meshsize  is  fine  enough, 
i.e.,  h  <  ho,  (assuming  that  k  is  related  to  ft  in  a  fixed  manner),  then  one  can  bound 


^''(n)  <  C{T} 


for  aU  nk  <  T, 


such  that  h  <  ho 


(3.65) 


this  bound  is  independent  of  the  mesh  size.  Thus, 

l|eS,^ll  <  7’’"+'  [C{T)r+^  $^+i(n)  (3.66) 

If  Bm+i  =  0(hP’"+i),  then  as  mesh  is  refined  the  error  in  the  LDS  approximation  decreases  like 

At  first  glance  the  discretization  (3.39)  is  a  first  order  in  time  explicit  scheme.  However,  any 
single  stage  explicit  or  implicit  discretization  of  any  order  can  be  brought  to  a  similar  form  once  the 
source  term  is  appropriately  redefined.  Therefore,  Theorem  2  can  be  modified  and  generalized 
to  yield  a  similar  result  for  a  general  discretization.  Such  a  generalization  would  be  futile  as 
finding  an  explicit  form  for  A,  B  and  B-m  might  require  matrix  inversion  or  computing  a  matrix 
polynomial.  Thus,  in  those  cases  the  bound  (3.52)  is  hard  to  compute. 

In  the  fuUy  discrete  setting,  as  in  the  semi-discrete,  an  inspection  of  Bm+i  might  indicate 
about  the  applicability  of  the  LDS  method  for  a  given  discretization. 


4  Stability,  consistency  and  convergence 

In  the  previous  section  a  method  for  obtaining  highly  accurate  approximations  using  an  expanded 
system  of  lower  order  approximations  was  analyzed.  It  should  be  proved  that  if  the  original  scheme 
was  consistent  and  stable  then  so  is  the  resulting  LDS  scheme. 

For  simplicity  of  presentation  the  discussion  is  limited  to  LDS  of  degree  one.  The  generalization 
to  LDS  of  general  degree  m  is  straightforward.  It  is  further  assumed  that  the  source  term  F  =  0. 
This  assumption  does  not  effect  stability  analysis  [14]  and  its  effect  on  consistency  will  be  shortly 
discussed. 
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In  an  LDS  of  degree  one  instead  of  using  the  stable  and  consistent  scheme, 


Un+l 

=  {I  +  kA)U^ 

(4.1) 

=  Uq 

(4.2) 

the  following  LDS  system  is  employed, 

KP  = 

{I+kA)Vo^^,  +  kV,^^, 

(4.3) 

KP  = 

(/  +  kA)Vr,, 

(4.4) 

Vo”.  = 

Uo 

(4.5) 

v.”.  = 

BiUq 

(4.6) 

where  =  B  is  defined  in  (3.39). 

The  stability  of  the  scheme  (4.3)-(4.4)  follows  from  the  structure  of  this  system.  It  consists  of 
a  principal  part  which  is  coupled  through  lower  order  terms.  The  principal  part  is  diagonal;  thus, 
its  discretization  is  stable  for  the  same  time  step  as  the  original  scheme.  The  lower  order  term 
does  not  aflFect  stability. 

The  consistency  of  the  single  equation  discretization  implies  that  the  LDS  is  a  consistent 
approximation  of  the  system 


Ut  =  Lu  +  t  (4.7) 

Tt  =  Lt 

for  initial  data 

u(a:,0)  =  uo(a;) 
r(a;,0)  =  to{x) 

In  the  LDS  approximation  =  BiUq.  Thus,  if  the  following  relation  holds 

Bi  =  0{hJ’^  +  ^^1)  with  P\,q\  >  1 

then  the  LDS  solution  is  a  consistent  approximation  of  (4.7)-(4.8)  with  initial  data 

u(x,0)  =  uo(2;) 
r(a;,0)  =  0 

Thus,  it  consistently  approximates  the  equation 

Ut  -  Lu  (4.11) 

u(a;,0)  =  Uo(®) 

The  source  term  for  the  equation  is  Bi  F.  Thus,  by  the  above  reasoning,  if  (4.9)  holds  then 
the  r  equation  is  a  consistent  approximation  of  an  equation  with  F  =  0.  Therefore,  this  term 
does  not  effect  consistency,  either. 

It  follows,  by  Lax  equivalence  theorem  [14],  that  the  consistency  and  stability  of  the  origi¬ 
nal  equation  together  with  condition  (4.9)  ensure  the  LDS  solution  convergence  to  the  analytic 
solution. 


(4.8) 

(4.9) 

(4.10) 
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The  LDS  approximation  maintains  the  stability  of  the  original  scheme,  however,  it  does  not 
necessarily  preserve  the  time-stability  of  the  discretization,  i.e.,  the  LDS  solution  might  exhibit  a 
non-physical  growth  in  time  although  the  underlying  discretization  did  not  allow  such  a  growth. 
This  phenomena  is  demonstrated  in  Lemma  1.  Recall,  if  is  a  discrete  operator  acting  on 
grid  functions  on  a  periodic  or  infinite  domain,  then  the  symbol  of  A^,  A^{d),  is  defined  by  the 
identity, 

^  ^h^ie-x  (4  J2) 

Lemma  1 :  Let  A^  be  a  space  discretization  of  a  scalar  time-dependent  equation  with  periodic 
boundary  conditions,  such  that  its  symbol  satisfies  A^{6o)  =  0,  for  7^  0  and  B^{9q)  7^  0.  Then 
the  error  in  the  LDS  system  based  on  A^  grows  polynomially  in  time  for  this  Fourier  component. 
The  order  of  the  polynomial  equals  the  LDS  degree. 

Proof:  Consider  a  semi-discrete  LDS  approximation  of  degree  one  based  on  A^.  According 
to  Section  3.1.1  it  has  the  form, 


u’}  =  A^u^  +  r'^  (4.13) 

=  A^ 

with  initial  data 

w^(a:,0)  =  Uq{x)  (4-14) 

T^{x,0)  —  B^Uq{x) 

The  solution  of  this  system  for  the  component  is, 

u^{0o,  t)  =  (u^(0o)  +  tB'^(9o)A^o))  =  (l  +  t  B'^(9o))  n\9o)  (4.15) 

The  proof  for  higher  degree  LDS  is  similar. 

Note,  that  the  property  A^(0o)  =  0  is  common  to  central  discretizations  of  the  first  derivative 
on  non-staggered  grids,  for  =  ’>'  (see  also  Sec.  6.3). 

It  should  be  emphasized  that  although  the  LDS  transformation  preserves  stability,  the  resulting 
algorithm  might  not  be  stable.  Consider,  for  instance,  a  discretization  satisfying  the  condition  of 
Lemma  1;  then  the  LDS  solution  for  the  component  wiU  grow  exponentially  with  the  number 
of  visits  to  the  fine  grid.  Other  possible  sources  for  such  a  growth  are  large  errors  introduced 
by  the  intergrid  transfer  employed  by  the  algorithm  which  are  not  damped  during  the  cycle  (see 
Sec.  5.1.4)  or  improper  use  of  Richardson  extrapolation  (see  Sec.  5.2).  These  last  remarks  can  be 
understood  once  the  LDS  algorithm  is  presented,  in  the  next  section. 

5  Large  Discretization  Step  (LDS)  Methods 

The  LDS  approximation,  introduced  in  the  previous  sections,  approximates  a  high  accuracy 
scheme  by  an  enlarged  system  of  equations  of  a  lower  accuracy  discretization.  In  the  present 
work  the  two  schemes  are  the  same  discretization  of  a  differential  operator  on  two  different  grids. 
An  LDS  algorithm  computes  the  fine  grid  solution  on  the  coarse  grid  by  solving  there  an  extended 
system  of  equations  which  are  initiaUzed  using  the  fine  grid.  The  LDS  system  is  integrated  on 
the  coarse  grid;  hence,  the  accuracy  of  the  correction  terms  deteriorates  at  a  rate  determined  by 
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that  grid  discretization.  However,  since  the  norm  of  these  terms  is  significantly  smaller  than  the 
solution  norm,  they  can  be  effectively  used  for  many  coarse  grid  time  steps.  Then,  the  fine  grid 
has  to  be  revisited  to  compute  new  initial  data  for  the  correction  terms.  These  terms  initialization 
is  computationally  costly  and  the  LDS  approximation  consists  of  more  equations  than  the  original 
problem;  however,  the  large  number  of  steps  performed  on  the  coarse  grid  before  revisiting  renders 
the  resulting  algorithm  very  efficient. 

This  section  details  the  algorithmic  implementation  of  these  ideas.  The  algorithm  is  presented 
both  in  its  Correction  Scheme  and  Full  Approximation  Scheme  forms  and  efficient  initialization 
procedures  for  these  schemes  are  described  and  analyzed.  Richardson  extrapolation  is  introduced 
to  the  LDS  method  which,  for  smooth  solutions,  yields  a  higher  order  approximation.  The  effi¬ 
ciency  of  the  algorithm  is  discussed  and  evaluated. 

Given  a  system  of  hyperbolic  differential  equations  with  coefficients  which  may  depend  on  x 
but  not  on  t,  of  the  form 

?£p^l-A(x,±)U(x,t)  =  F(x) 

Mu{x,t)  =  0 

17(a;,0)  =  Uo{x) 

where  C  iR^  and  ^  =  (gl^,  gf^). 

Consider  a  discretization  of  the  form 

=  E(x,k,h)U^  +  S{x,k,h)F^  (.5.2) 

where  h,k  denotes  Aa;  and  At,  respectively,  and  17"  =  Uf  approximates  U{jh,nk).  In  this  work 
E{x,k,h)  is  an  explicit  or  implicit  two  level  time  marching  operator.  However,  the  method  may 
be  used  with  multilevel  integration  schemes,  as  weU.  In  the  sequel,  the  notation  E^'’^  will  be  used, 
omitting  the  possible  dependence  on  x.  In  an  LDS  application  two  grids  (in  space- time)  are  given, 
a  fine  one  with  spacing  {h,  k)  and  a  coarse  with  (if.  A'),  where  H  =  ah,  K  =  ak.  Given  (x,  0) 
on  the  {H,  K)  grid,  one  needs  to  calculate  the  solution  up  to  a  prescribed  final  time  T  and  obtain 
the  fine  grid  accuracy. 

5.1  The  LDS  Method  of  General  Degree 

The  degree  of  an  LDS  approximation  is  defined  as  the  number  of  its  correction  terms.  The 
error  bounds  obtained  in  Section  3  suggest  that  in  many  cases  (e.g.,  for  smooth  solutions)  the 
higher  the  degree  the  better  the  LDS  approximates  the  fine  grid  solution.  In  practice,  efficiency 
considerations  limit  the  degree  to  at  most  two  (see  Section  5.3.2).  In  the  sequel,  algorithms  for 
LDS  of  general  degree  are  described. 

The  initialization  procedures  necessitate  the  transfer  of  the  solution  from  the  coarse  grid  to 
the  fine  and  back.  In  this  section,  it  is  assumed  that  appropriate  intergrid  transfers  are  given. 
The  order  and  properties  of  these  transfers  will  be  discussed  in  Section  6.1.  For  presentation 
simplicity  the  algorithm  is  described  for  the  case  F{x)  -  0;  the  treatment  of  a  source  term  is 
straightforward. 


for  X  €  D,  t  e  [o,r] 

for  X  £  do.  (5.1) 

for  X  e 
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5.1.1  Correction  Scheme  LDS 

The  fundamental  idea  of  the  LDS  method  is  to  look  for  correction  terms  to  the  coarse  grid 
ei^uations,  derive  and  solve  the  ec^uations  satisfied  by  these  terms  to  obtain  the  fine  grid  accuracy 
on  the  coarse  grid.  In  Section  3,  it  was  shown  that  for  linear  problems  these  terms  approximately 
satisfy  the  same  equation  as  the  solution.  The  resulting  system  of  equations,  which  is  valid  only 
for  linear  problems,  was  named  Correction  Scheme  LDS. 

Correction  Scheme  LDS  of  Degree  One  when  H  =  2h.  For  clarity  of  presentation,  the 

algorithm  is  first  presented  its  simplest  form,  i.e..  Correction  Scheme  LDS  of  degree  one  with 

H  _  ^  _  o 

h  -  k  ^  ^ 

The  algorithm  consists  two  stages  :  initialization  of  the  correction  terms  using  the  fine  grid, 
and  time  marching  on  the  coarse  grid  for  a  predetermined  number  of  steps.  The  results  in  Section  3 
show  that  for  linear  problems  the  correction  terms  satisfy  approximately  the  same  equation  as  the 
solution.  However,  they  do  no  indicate  how  to  effectively  and  efficiently  compute  initial  values  for 
these  terms.  The  requirement  that  the  LDS  solution  should  yield  the  fine  grid  solution  suggests 
that  on  the  first  time  steps  these  solutions  should  be  identical.  This  observation  leads  to  the 
following  initialization  routine. 

Initialized^^, 


Solve  tf^+f  =  uN+^  ^ 


II 

_  jH,K  jjN+l 

h,k 

-  kS*' 

VoT 

= 

+ 

N 

=  iV  + 

1 

The  initialization  consists  of  interpolating  the  LDS  solution  to  the  fine  grid  and  stepping  for 
the  same  time  on  the  fine  and  coarse  grids.  Then,  the  correction  term  is  set  to  the  discrepancy 
between  the  two  solutions  and  the  LDS  solution  is  updated  to  the  fine  grid.  Here,  stands 

for  an  intermediate  value  assigned  to  this  variable. 

The  time  advance  of  the  LDS  has  the  following  simple  form, 

LDS  Method  of  Degree  One,  Correction  Scheme 

Initialize 
iV  =  0 

While  N  <  f  J]  Do 

Call  Initialize(Fo^i,Fi^) 

For  i  =1,. . .  ,Revisit  Do 
Solve 
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Set  jV  =  iV  +  1 

End 

End 

The  relation  between  the  theory  derived  in  Section  3  and  the  initialization  procedure  is  ana¬ 
lyzed  in  Section  5.1.2. 

Correction  Scheme  LDS  of  General  Degree  when  H  =  2h.  The  correction  term  in  the 
LDS  of  degree  one  is  initialized  to  equalize  the  fine  and  LDS  solutions.  If  one  could  solve  the 
exact  equation  satisfied  by  this  term,  which  is  approximately  the  fine  grid  equation,  then  the  fine 
grid  and  the  LDS  solutions  would  be  identical.  Clearly,  this  would  be  as  difficult  as  obtaining 
the  fine  grid  solution  on  the  coarse  grid.  Instead,  one  can  add  a  new  term  to  correct  the  first 
correction  term  equation.  Thus,  in  an  LDS  of  general  degree,  the  i  correction  term  may  be 
viewed  as  correcting  the  {i  —  1)*^  term. 

The  general  degree  algorithm  consists  of  two  stages  :  initialization  of  the  correction  terms 
and  time  marching  on  the  coarse  grid  for  a  predetermined  number  of  steps.  The  time  marching 
procedure  has  the  following  simple  form,  where  term  i  corrects  term  f  -  1  for  1  <  *  <  d. 

LDS  Method  of  General  Degree  -  d  ,  Correction  Scheme 

Initialize 

N  =  0 

While  N  <  [^1  Do 

Call  InitializeCFo^, ...,'Cj^,d  ) 

For  i  =  1,. .. .Revisit  Do 

Solve 

For  1  =  d-l,...,0,  Step  =  -1,  Do 

sol,.  VS  +  VSti 

End 

Set  N  ^  N  +  1 

End 

End 

In  an  LDS  of  degree  d,  the  correction  term  (1  <  i  <  d)  corrects  the  (i  - 1)*'"  term  and  both 
satisfy  approximately  the  same  equation  as  the  coarse  grid  solution.  Therefore,  the  initialization 
of  the  term  to  correct  the  (z  —  term  is  identical  to  the  initialization  of  the  first  term  to 
correct  the  solution  in  the  LDS  of  degree  one.  Once  is  initialized  and  is  updated  using 

this  value;  the  lower  index  variables  can  be  time  advanced  using  these  new  values.  This  procedure 
is  repeated  for  all  1  <  i  <  d.  Section  5.1.2  outlines  the  connection  between  the  error  bound 
derived  in  Section  3.1.1  and  these  intuitive  arguments  which  are  implemented  in  the  procedure 
listed  below. 

Initialize(F(j^,  ...,Fj^,d  ) 
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wf  =  Ehu^  + 


Figure  1:  Initialization  of  an  LDS  of  degree  two.  With  the  notation  :  =  Vo,2,  =  Ii,2! 

=  ^2,2 
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Correction  Scheme  LDS  of  General  Degree  when  H  =  ah.  The  above  described  proce¬ 
dure  can  be  easily  adapted  for  the  general  case  when  ^-  =  a.  The  simplest  approach  is  to  perform 
a  time  steps  on  the  fine  grid  for  each  coarse  grid  time  step  and  initialize  the  correction  terms 
correspondingly,  see  Figure  2.  This  procedure  of  direct  initialization  is  very  costly  and  greatly 


Figure  2:  Direct  Initialization 

reduces  the  LDS  efficiency. 

In  case  a  is  a  composite  number,  e.g.,  a  =  2^,  a  more  efficient  approach  is  available,  exploit¬ 
ing  the  LDS  high  accuracy  by  employing  intermediate  grids.  In  this  approach,the  simultaneous 
initialization,  the  fine  grid  is  used  to  initialize  an  LDS  system  of  degree  one  on  the  grid  Hi  =2h] 
since  this  approximation  is  very  accurate  it  can  be  used  to  initialize  an  LDS  of  degree  one  on  grid 
H2  =  2 Hi-  This  procedure  is  repeated  until  the  correction  term  on  grid  Hi  =  H  is  initialized. 
This  process  is  repeatedly  employed  for  all  correction  terms.  In  this  method,  a  correction  term 
on  a  coarse  grid  is  initialized  as  soon  as  enough  time  marching  on  finer  grids  was  performed,  see 
Figure  3. 

A  few  important  points  should  be  emphasized  regarding  the  initialization  procedures.  First, 
in  this  work  the  computationally  efficient  method  was  favored.  There  is,  however,  a  trade  off 
between  the  computational  cost  and  the  storage  requirements  of  these  two  methods  (see  Sec.  -5.3). 
Therefore,  whenever  storage  is  limited,  direct  initialization  might  be  preferred.  Second,  it  might 
have  been  expected  that  direct  initialization  will  yield  more  accurate  solutions  than  the  simulta¬ 
neous  approach.  However,  our  numerical  results  show  that  these  methods  are  indistinguishable 
for  integrations  times  of  interest,  i.e.,  as  long  as  the  error  in  the  LDS  solution  is  small.  Last, 
the  direct  initialization  is  of  practical  interest,  being  used  with  Richardson  extrapolation  (see 
Sec.  5.2). 

A  simple  way  to  predict  the  LDS  performance  is  to  look  at  the  relative  magnitude  of  the 
correction  terms  immediately  after  initialization.  According  to  the  result  presented  in  Section  2, 
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Figure  3:  Simultaneous  Initialization 


the  ratio  should  be  roughly  constant.  Thus,  a  large  variation  in  this  quantity  suggests 

a  large  error  in  the  initialization  of  Tj,  causing  the  LDS  failure. 

The  error  bounds  derived  in  Section  3  apply  to  approximations  of  general  degree  m.  In  practice 
however,  due  to  efficiency  considerations  two  is  the  highest  degree  used  (see  Sec.  5.3.2). 

5.1.2  Initialization  Analysis 

The  procedure  for  the  correction  terms  initialization  was  justified  by  the  intuitive  arguments  that 
when  properly  initialized,  the  LDS  solution  should  agree  with  the  fine  grid  solution  for  the  first 
time  steps;  and  that  in  an  approximation  of  degree  d  the  term  corrects  the  (i  -  1)*^ 
term  (for  1  <  i  <  d).  Thus,  each  term  is  initialized  similarly  to  the  manner  the  first 

term  is  set  to  correct  the  solution.  This  section  outlines  the  relation  between  the  error  bound 
derived  in  Section  3.1.1  for  the  semi-discrete  case  with  the  commutativity  assumption  and  the  way 
initialization  is  implemented.  Specifically,  it  will  be  shown  that  the  solution  at  the  termination  of 
the  initialization  procedure  described  in  Section  5.1.1  and  the  solution  of  the  LDS  approximation 
introduced  in  Section  3.1.1  are  equal,  up  to  higher  order  terms  and  multiplicative  constants  which 
are  used  for  computational  efficiency. 

For  simplicity,  the  analysis  is  performed  for  the  case  y  =  2.  It  is  assumed  that  the  fine 
and  coarse  grid  spatial  discretizations  commute,  i.e.,  [L^,  L^]  =  0;  and  that  intergrid  transfers 
introduce  no  error.  By  an  abuse  of  notation  these  transfers  are  omitted  from  the  analysis  in  this 
section.  Nevertheless,  whenever  the  fine  and  coarse  grid  solutions  appear  in  the  same  formula,  it 
should  be  understood  that  an  appropriate  restriction  of  the  fine  grid  solution  to  the  coarse  grid 
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is  implied.  Similarly,  whenever  the  fine  and  coarse  grid  operators  appear  in  the  same  formula 
it  stands  for  a  restriction  of  the  fine  grid  operator  to  the  coarse  grid.  The  properties  of  these 
transfers  necessary  to  guarantee  the  algorithm  performance  are  discussed  .in  Section  6.1. 

In  Section  3.1.1  an  error  bound  was  derived  for  a  semi-discrete  LDS  approximation  of  general 
degree  d.  This  approximation  which  is  given  by  equations  (3.3)-(3.4)  can  be  succinctly  written 


L^'Vi,d  +  Vi+i^d  foT  i  <0  <  d 


=  Vd,d 


with  initial  conditions 


Vi^i{x,0)  —  uo{x) 


where  uo{x)  is  evaluated  at  grid  points.  The  solution  of  this  system  is  given  by, 


for  0  <  i  <  c? 


For  an  LDS  algorithm  of  degree  one  denote  by  the  variables  approximating 

respectively.  At  time  At  the  solution  of  the  LDS  approximation  is  given  by, 

vo,i{x,At)  =  [l  +  At[L>^-L^)]e^''^^uo{x)  (5.6) 

At)  =  uo(a:) 

For  linear  problems,  the  time  marching  operator  approximates  Therefore,  in  the  fol¬ 
lowing  semi-discrete  analysis  will  denote  the  time-stepping  on  grid  H  for  time  At.  Assuming 

u{x,  At)  =  0(1),  then,  at  the  end  of  the  initialization  phase  the  variables  satisfy 

«ii)s(a;,At)  =  e^^^^uo{x)=  I -\- At  {l^  -  -  L^Y  e^^^^uo(a;) -|- h.o.t 

=  ^0,1(2^)  ^0  +  0  {{Atf  -  L^Y')  (5-'^) 


-(a;,At)  = 


gL'-At  _  gL»Ai 


)  «o(a;) 


=  At  I  +  ^  -  l^Y]  e^^'^^uoix)  -f  h.o.t 

=  At  |ui,i(a;,At) -f  0  ^At  (x*  - 


The  factor  At  is  maintained  to  reduce  the  number  of  multiplications  during  the  time-stepping 
stage.  It  can  be  seen  that,  up  to  higher  order  expressions,  the  algorithm  correctly  initializes  the 
first  correction  term. 


22 


For  an  LDS  algorithm  of  degree  two  denote  by  ^  variables  approximating  vo,2,  ^1,2 

and  V2,2-,  respectively.  At  time  2At  the  solution  of  the  LDS  approximation  is  given  by, 

vo,2{x,  2At)  -  1  +  2At  +  2{Aif  j 

ui,2(x,2At)  =  [l  +  2At{L^-L^)\e^‘'^^^[L^-L^)uo{x)  (5.9) 

r;2,2(a:,2At)  =  uq{x) 

It  can  be  easily  seen  that  at  the  end  of  the  initialization  procedure  the  variables  r,  (p  satisfy, 

u^^,{x,2At)  =  e^''^^^UQ{x)  =  v^,2{x,2At)  +  o({Atf[L^-L^Y^  (.5.10) 

r(2:,  2At)  =  uo{x)  = 

=  At  {(l'^  -  ^  +  h.o.t 

=  At|ui,2(a:,2At)  +  0  (^At(x^-L^)^^l  (.5.11) 

If  {x,  2 At)  = 

=  {Atf  I  -L«)\  At  {l^  -  L^y  +  L^)  }  e^^^^\o{x)  +  h.o.t. 

=  (A/)2  |u2,2(x,  2A/)  +  0  (At  (x'^  -  X^)^) I  (5.12) 

It  can  be  seen  the  algorithm  correctly  initializes  the  correction  terms  in  the  LDS  of  degree  two, 
up  to  higher  order  terms  and  multiplicative  constants. 

5.1.3  Full  Approximation  Scheme  LDS 

The  FAS  form  of  the  LDS  is  appropriate  for  both  linear  and  nonlinear  problems.  However,  since 
the  Correction  Scheme  has  a  simpler  form  and  necessitates  less  modifications  to  the  code,  it  is 
more  conveniently  used  for  linear  problems. 

In  the  present  work  only  LDS  in  FAS  form  of  degree  one  was  implemented  ,  and  it  will  be 
described  for  a  homogeneous  system  of  equations. 

Recall,  the  Full  Approximation  Scheme  of  degree  one  is  given  by, 

u'l  =  p”{u'^)  +  v^  -u^  (5.13) 

where  may  be  either  a  linear  or  nonlinear  operator  and  +  r,  with  r  corresponding  to 

the  first  correction  term  in  the  Correction  Scheme  form. 

The  algorithmic  implementation  is  slightly  more  involved  than  the  Correction  Scheme  as  it 

^  — TT  TC 

requires  some  modification  of  the  time  marching  procedure.  Denote  by  ’  the  integration 
scheme  obtained  by  modifying  the  coarse  grid  operator  to  time  advance  the  LDS  system  of 
degree  one  (5.13). 
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The  algorithm  has  the  following  form 

IDS  Method  of  Degree  one  ,  Full  Approximation  Scheme 


Initialize 

N  =  0 

While  iV  <  [f]  Do 

Call  Initialize(Fo4,Ti^i) 

For  i  =1,. . .  jRevisit  Do 
Solve  Vf , 

Set  iV  =  JV  +  1 

End 

End 

Here,  Fj  denotes  the  vector  (Fo,i,  Fi^i).  For  presentation  simplicity,  the  initialization  procedure 
is  described  for  the  case  x  =  2.  Generalizations  to  the  case  f-  =  2^  are  identical  to  those  for  the 
Correction  Scheme.  In  Section  5.1.2  it  was  shown  (see  Eq.  (5.8)-(5.7)),  that  after  the  Correction 
Scheme  initialization  contains  the  fine  grid  solution  and  F/^  contains  Atr(a;,  At). 

Thus,  the  FAS  initialization  consists  of  the  Correction  Scheme  initialization  supplemented  with 
the  additional  computation  of  +  r)(a;,  At)  at  the  end  qf  the  procedure.  For  completeness,  the 
whole  initialization  procedure  for  the  Full  Approximation  Scheme  is  listed  below 

Initialized^,  Fj^) 


Set 

Solve  uN+i2=i  ^ 

F(J+i  =  E^’^  Fo^ 

Set  Fi^i+^  =  _  yN+i 


Set 


N 


=  A  +  1 


m=l  ,2 


The  generalization  to  higher  degree  LDS  is  straightforward,  instead  of  the  original  variables 
(Fo,rf,  Fi,d, . . . ,  Vd,d)  a  new  set  of  variables  (Vo,d,  Vo,d  +  Vi,d,  ■  ■  ■,Vo4  +  - ■  ■  +  Vd^d)  is  introduced.  The 
equations  satisfied  by  these  new  variables  can  be  determined  from  the  equivalence  between  the 
Full  Approximation  Scheme  and  the  Correction  Scheme  for  linear  problems. 

5.1.4  Treatment  of  Boundary  conditions 

The  LDS  treatment  of  the  boundary  conditions  will  be  discussed  under  a  restrictive  commutativity 
assumption,  which  at  this  stage  we  do  not  know  how  to  dispose. 
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(5.14) 


Assume  satisfies  the  boundary  condition 

Mu^  +  Nu^  =  g 

If  [M,  iV]  =  0,  then  t  =  Nu^  satisfies, 

Mt  =  M  N  =  N  Mvf^  =  Ng  -  N'^u  (5.15) 

If  Mr  <  r  the  right  hand  side  term  may  be  neglected  and  r  satisfies  the  boundary  condition 

Mr  =  0  (5.16) 

Otherwise,  when  higher  degree  LDS  is  employed,  the  rj  satisfy  the  boundary  conditions 

Mri  +  r2  =  fifi  (5.17) 

MTk  =  0  (5.18) 

where 

Tj  —  N^t  (5.19) 

(5.20) 

Thus,  a  large  error  at  the  boundary  discretization  may  require  adding  correction  terms  and 
corresponding  equations  in  the  whole  domain. 

The  assignment  of  the  appropriate  boundary  conditions  to  the  correction  terms  when  the 
commutativity  assumption  does  not  hold  should  be  further  investigated. 

The  presence  of  non-periodic  boundary  conditions  may  pose  problems  even  when  the  commu¬ 
tativity  assumption  holds.  This  is  due  to  errors  introduced  by  the  one-sided  high  order  interpola¬ 
tion  near  the  boundaries.  These  large  and  localized  errors  excited  during  the  initialization  phase 
might  not  be  damped  before  the  next  visit  to  the  fine  grid;  resulting  in  an  error  which  grows 
exponentially  in  the  number  of  visits  to  the  fine  grid.  At  this  stage  of  research,  it  seems  that 
one  should  use  the  differential  equation  to  design  appropriate  near  boundary  interpolation  with 
reduced  errors  (see  Sec.  7  for  an  example). 

5.2  Richardson  Extrapolation 

The  simultaneous  time  stepping  on  two  grids  during  the  initialization  phase  can  be  used  to 
estimate  the  local  truncation  error.  This  estimate  can  be  used  in  various  ways.  In  [1,  2,  3]  it  was 
used  to  implement  adaptive  mesh  refinement  for  hyperbolic  equations.  In  the  multigrid  method 
this  estimate  is  used  for  the  r  extrapolation  technique  which  applies  a  weighted  transfer  of  the 
correction  term  to  obtain  higher  accuracy  using  a  lower  order  scheme  [4]  . 

For  simplicity,  let  Qh  be  a  two-level  explicit  difference  operator.  If  the  solution  is  smooth 
enough,  the  local  truncation  error  is 

u{x,t  +  k)  —  Qhu{x,t)  =  k[k'^^a{x,t)  +  h'^‘^b{x,t)]  + (5.21) 

=  T  +  kO{k’^^+'^  +  (5.22) 
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where  the  leading  term  is  denoted  by  r.  If  tt  is  smooth  enough,  then  if  one  takes  two  time  steps 
with  the  method  Qh-,  the  leading  error  is  2r.  That  is, 

u{x,  t  +  2k)  -  QIu{x,  t)  =  2t  +  (5.23) 

Let  Q2h  be  the  same  difference  scheme  as  Qh  but  based  on  mesh  width  of  2h  and  2k.  Also,  assume 
that  the  order  of  accuracy  in  time  and  space  are  equal,  qi  =  q2  =  q.  Then 

u{x,t  +  2k)  -  Q2huix,t)  =  {2k)[{2kya{x,t)  + {2h)^b{x,t)]  + (5.24) 

=  2«+V  +  0(/i»+2)  (525) 

Since  u{x,t  +  2k)  —  Q^u{x,t)  ss  2r,  forming  the  difference 

gives  an  estimate  of  the  local  truncation  error  at  time  t.  In  other  words,  the  difference  between 
the  solution  on  grid  {2h,  2k)  and  (h,  k)  uses  to  estimate  of  the  local  truncation  error. 

This  procedure  has  several  advantages.  First,  it  is  not  necessary  to  know  the  exact  form  of 
the  truncation  error  to  apply  it.  The  error  estimation  procedure  is  independent  of  the  difference 
method.  The  restriction  of  this  method,  that  the  accuracy  in  time  and  space  should  be  the  same, 
is  not  a  severe  one.  Many  popular  finite  difference  methods  share  this  property,  for  example, 
second  order  methods  like  Lax  Wendroff  or  MacCormack’s  method  and  Leap  Frog,  and  first  order 
method  such  as  upstream  differencing.  For  methods  where  the  accuracy  in  space  and  time  is  not 
the  same,  a  more  expensive  variant  of  this  procedure  is  possible.  For  example,  one  could  estimate 
the  spatial  and  temporal  error  separately:  first  keep  k  constant  and  take  a  step  based  on  2h 
differences,  then  keep  h  constant  and  take  a  step  with  time  step  2k.  Other  variations  are  possible. 
In  the  present  work  the  LDS  with  Richardson  extrapolation  was  employed  only  for  schemes  with 
the  same  spatial  and  temporal  accuracy.  The  usefulness  of  this  approach  applied  to  discretizations 
without  this  property  should  be  further  investigated. 

The  initialization  step  of  the  LDS  method  computes  the  term  Q\u  -  Q2hU  and  uses  it  as  the 
initial  value  of  the  correction  equation  on  the  next  coarser  grid.  Taking 

{QU  -  Q2hu)  (5.27) 

as  the  initial  value  will  yield  an  (^(/i®"*"^)  approximation  on  the  coarse  grid. 

The  initialization  of  the  extrapolated  LDS  of  degree  one  is  performed  directly  from  the  finest 
grid.  The  extra  cost  associated  with  this  initialization  is  compensated  by  the  added  accuracy. 

Richardson  extrapolation  is  based  on  Taylor  expansion  of  the  error  and  is  valid  only  for  smooth 
data.  For  non-smooth  solutions,  this  procedure  is  incorrect  and  might  lead  to  an  error  that  grows 
exponentially  with  the  number  of  visits  to  the  fine  grid.  Therefore,  a  great  care  should  be  taken 
when  considering  this  method,  to  ensure  that  the  discretization  is  dissipative  enough  to  prohibit 
any  undesired  growth.  Nevertheless,  since  the  dissipative  schemes  employed  by  the  LDS  damp  the 
oscillatory  components  (see  Sec.  6.2),  this  technique  might  give  exceffent  results  when  properly 
used  (see  Figure  12). 
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5.3  Work  Considerations 

The  amount  of  computational  work  and  the  memory  requirement  in  a  cycle  of  an  LDS  of  degree  m 
will  be  evaluated  and  compared  with  the  corresponding  requisites  on  the  ;finest  grid  in  the  cycle. 

The  simultaneous  time-stepping  on  the  finer  grids  during  initialization  constitutes  a  large 
fraction  of  the  algorithm  computational  cost  and  dominates  its  storage  requirements.  There  is  a 
trade  off  between  the  storage  requirements  and  the  efficiency  in  the  two  initialization  procedures 
described  in  Section  5.1.  In  the  present  work  the  computationally  efficient  simultaneous  initial¬ 
ization  was  employed  since  only  moderate  storage  was  needed.  However,  when  storage  is  limited, 
direct  initialization  might  be  favored.  In  this  section,  only  the  efficiency  of  the  simultaneous 
scheme  is  analyzed,  but  the  storage  requirements  are  compared  for  both  methods. 

In  order  to  simplify  analysis,  it  will  be  assumed  that  the  finest  and  coarsest  grid  meshsizes 
satisfy  H  =  2^h;  and  that  a  grid  is  refined  by  halving  its  meshsize. 

The  problem  is  solved  in  a  d-dimensional  space,  for  d  =  2, 3.  Typically,  real  world  problems 
occur  in  3-dimensional  space. 

5.3.1  Storage  Requirements 

The  LDS  method  of  degree  m  employs  on  the  coarsest  level  m  +  1  times  as  many  equations  as 
on  the  finest  grid;  while  on  intermediate  grids,  twice  the  number  of  the  fine  grid  equations  are 
solved.  The  number  of  points  in  a  spatial  grid  on  any  level  is  2“^  larger  than  in  the  next  coarser 
one.  In  the  simultaneous  initialization,  see  Figure  3,  if  all  grids  are  allocated  simultaneously  the 
storage  requirement  is, 

2^ j  ^ 

where  S  is  the  storage  required  for  the  finest  grid.  Noting  that  at  all  times  merely  two  grids  are 
time  advanced  simultaneously,  then  a  careful  management  of  memory  may  reduce  this  requirement 
to 

^  ^  (^-29) 

These  two  requisites  are  equivalent  for  /  <  2  (i.e.,  for  ^  =  2,  or  4)  as  is  the  case  in  the  present 
research.  Therefore,  this  possible  small  reduction  of  memory  usage  will  not  be  further  elaborated. 

If  memory  is  at  premium,  storage  may  be  traded  for  efficiency  by  using  direct  initialization^ 
see  Figure  2.  This  procedure  employs  only  the  finest  and  coarsest  grids  with  memory  requirement 
of 

(1  +  ^)^  (5-30) 

It  should  be  noted  that  typically  I  =  2,  d  >  2,  and  due  to  efficiency  considerations  the  degree 
satisfies  m  <  3  (see  Section  5.3.2);  thus,  the  storage  overhead  associated  with  the  LDS  algorithm 
is  fairly  small. 

5.3.2  Efficiency 

The  computational  cost  of  fine  grid  time-stepping  relative  to  the  cost  of  obtaining  the  same 
solution  using  the  LDS  algorithm  will  be  evaluated.  In  this  estimate,  the  cost  of  the  intergrid 
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transfers  is  neglected,  as  in  many  cases  it  is  small  relative  a  fine  grid  time  step.  First,  the  work 
associated  with  the  initialization  of  an  LDS  of  degree  m  is  computed.  Define 

j  _  I  The  cost  to  initialize  an  LDS  of  degree  m  using  the  simultaneous 
I  initialization  when  the  spatial  dimension  is  d  and  H  =  2^h 

An  inspection  of  Figure  .3  leads  to  the  following  formula  for  Ii^m{d)-, 


T  /  ^  m(m  +  1)  4-^1 

Il,m{d)  =  2m+  (5..32) 

During  initialization  a  time  equal  to  m  2‘  fine  grid  time  steps  is  marched.  Denote  by  N  the 
number  of  coarse  grid  time  steps  performed  before  revisiting  the  fine  grid.  The  efficiency  of  an 
LDS  cycle  is  defined  as  the  computational  cost  of  time-stepping  on  the  fine  grid  relative  to  the 
work  required  to  obtain  the  same  solution  on  the  coarse  grid  with  an  LDS  cycle.  It  is  given  by 
the  formula, 

/,,„(d)  +  W(m+l)2-'“ 

The  following  table  and  figure  list  the  LDS  efficiency  for  H  =  ih  and  d  =  2, 3,  for  various  values 
of  A'. 


Revisit 

Efficiency  for  2D  problems 

Efficiency  for  3D  problems 

m  =  1 

m  =  2 

m  =  3 

m  =  1 

m  =  2 

m  =  S 

5 

6.98 

4.23 

4.25 

10 

10.83 

6.35 

Kin 

15 

13.15 

8.00 

5.82 

22.38 

12.36 

8.83 

20 

15.81 

9.32 

6.75 

25 

17.52 

10.41 

7.53 

32.80 

Typically  to  multilevel  methods,  this  algorithm  efficiency  increases  with  the  problem  dimen- 
sionality. 


6  Fourier  Analysis 

Fourier  analysis  is  a  major  tool  for  the  analysis  of  numerical  approximations  to  hyperbolic  equa¬ 
tions  [17],  as  well  as  for  analyzing  multigrid  algorithms  [4].  In  this  section  it  is  employed  to  obtain 
necessary  conditions  for  algorithm  convergence. 

6.1  Properties  of  the  intergrid  transfers 

The  transfer  of  the  solution  between  the  various  grids  plays  a  central  role  in  the  LDS  algorithm. 
In  the  initialization  stage  the  solution  is  first  interpolated  to  the  finer  grids  and  after  several  time 
steps  on  these  grids  is  restricted  back  to  the  coarsest  grid.  For  an  LDS  of  higher  degree,  this  pro¬ 
cedure  is  repeated  for  the  correction  terms.  Inevitably,  this  process  introduces  errors.  Therefore, 
an  appropriate  choice  of  these  operators  is  essential  to  guarantee  the  algorithm  performances  and 
the  desired  accuracy  of  the  solution. 
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Figure  4:  LDS  efficiency  for  2D  and  3D  problems  when  H  =  4h.  (a)  IDS  of  degree  one  (b)  LDS 
of  degree  two  (c)  LDS  of  degree  three. 

The  analysis  will  be  first  performed  for  first  order  operators.  This  is  no  limitation  since  every 
problem  may  transformed  to  a  first  order  system  by  introducing  additional  variables.  Assume 
that  ^  is  fixed  on  aU  grids;  thus,  the  temporal  error  can  be  expressed  in  terms  of  h  = 
Ax.  For  simplicity,  the  analysis  wiU  be  performed  for  one  dimensional  problems  with  H  =  2h; 
generalizations  to  higher  dimensions  is  straightforward.  Assume  that  the  spatial  discretization  is 
of  order  p.  Furthermore,  assume  that  the  order  of  the  spatial  discretization  coincides  with  the 
order  of  the  fuU  discretization.  Thus,  for  this  analysis,  the  effect  of  the  temporal  error  on  the 
accuracy  may  be  ignored  by  investigating  a  semi-discrete  system. 

Let  the  Ijf,  iff  be  the  interpolation  and  restriction  operators,  respectively;  and  let  1^(0),  lff{B) 
be  their  corresponding  symbols.  Assume  that  for  smooth  components  the  intergrid  operators 
satisfy, 

ifj{9)  =  1  +  c^Pi  (6.1) 

iff  {9)  =  1  +  C0P2 

where  throughout  this  section  c  is  a  generic  constant.  For  the  harmonic  oscillatory  component 
9'  =  9  4- holds, 

i'}j{9')  =  -1-0(0''*+^)  (6.2) 

iff  {9')  =  + <9(092+1) 

Denote  by  i{9)  =  iff{9)ijj{9).  Then  under  the  previous  assumptions,  for  smooth  components 

i{9)  «  {l  +  c9P^){l  +  c9P^)=  l  +  c(^Pi  -l-^P2)  +  h.o.t  (6.3) 
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For  the  harmonic  oscillatory  components  holds, 


7(0')  =  +  h.o.t  (6.4) 

Let  L^{6)  denote  the  symbol  of  L^.  For  LDS  of  first  degree  one  wants  to  compute, 

f{0,  At)  =  -  eL"(20)AiJ 

Instead  one  computes, 

+  ■ao(0)  +  h.o.t  (6.6) 

For  the  low  frequencies,  if  L  is  the  first  derivative, 

L\0)  =  i-  +  ch^{^-j  +h.o.t  (6.7) 

For  the  high  frequencies, 

L\e)\  =  0  (i)  (6.8) 

Thus,  for  this  operator, 

f(6',At)  =  c0P‘''^^uo(6>)  +  h.o.t  (6.9) 

In  order  to  obtain  the  desired  accuracy,  the  following  inequalities  should  hold 


ffpl 

< 

ap+1  ^ 

h 

QV2 

< 

h 

^91+92 

< 

h 

(6.10) 

(6.11) 

(6.12) 


In  general,  for  first  order  hyperbolic  equations,  ^  =  0(1),  thus  this  term  may  be  neglected  and 
one  obtains  conditions  on  the  order  of  the  intergrid  transfers  required  to  guarantee  the  algorithm 
performances 


P+2  <  pi,P2  (6.13) 

P  +  2  <  9i  +  ?2 

It  follows  from  the  above  bounds  that  if  ^  ^  0,  increasingly  higher  order  interpolations 
will  be  necessary.  Hence,  decreasing  the  time  step  on  a  fixed  spatial  grid  will  have  undesired 
consequences.  In  practice,  one  tries  to  use  a  time  step  as  large  as  possible,  thus,  this  observation 
poses  no  real  restrictions. 

If  L  is  an  operator  of  order  m,  the  previous  argument  implies  that, 

f((9.  At)  =  c  +  h.o.t  (6.14) 
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Therefore,  jn  order  to  obtain  the  desired  accuracy,  the  following  inequality  should  hold 


(6.15) 


and  simUar  inequalities  should  hold  for  p,,  (q,  +  q,).  This  implies  the  Mowing  conditions, 

p  +  m  +  l  < 

P  ^  "hi  ^  “h  ^2 

provided  ^  =  0(1),  i.e.,  it  is  bounded  away  from  zero. 

coJrt  r  “h"'  ‘™’  i  “  “““““ "“y  ^  “  w'Mscd  to 

correct  u  Hence,  tlmse  conditions  are  sufficient  for  IDS  of  degree  two.  as  well.  It  follows  that 
these  conditions  are  sufficient  for  an  LDS  of  a  general  degree  m. 

These  orders  of  the  intergrid  transfer  are  necessary  to  maintain  the  scheme  accuracy.  In  prac¬ 
tice,  one  rnight  prefer  to  employ  transfers  of  order  higher  than  the  minimum  necessary.  Consffier 
vhi^rTlf  ^  dissipative  scheme  which  strongly  damps  the  oscillatory  components.  Thus,  each 
o  e  fine  grid  somewhat  damps  the  smoother  components  of  the  solution  through  the  inter- 
polation  error  which  transfers  some  of  their  energy  to  the  fine  grid  oscillatory  components.  High 

order  interpolations  transfer  less  energy  to  these  components  and  therefore  might  be  preferable 
m  such  circumstances. 

Another  undesircd  property  of  the  intergrid  transfer  is  captured  in  the  following  lemma. 

+  'f'  1.  initialuationpr^u,^ 

introduces  an  0(1)  error  in  the  Fourier  component  0. 

Proof  :  By  the  assumption,  instead  of  f(0.  At)  in  Eq  (6.9)  one  computes 
[l(0)e^  ~  «o(^)  =  f{0,  At)  +  [c  (0^^  +  0f^  )e^'‘(^)^' 

+  (C  -  1  -  ^ 

The  staggered  grid  transfers,  often,  satisfy  the  condition  of  Lemma  2  for  the  osciUatory  com- 

haThtSe  oraT?  ""  discretizations  are,  usuaUy,  dissipative  (see  Sec.  6.2),  this'r 

nas  little  practical  consequences. 

6.2  The  Symbol  of  the  LDS  cycle 

The  LDS  method  employs  several  operators  which  should  be  simultaneously  analyzed  in  order  to 
ensure  the  proper  performance  of  the  algorithm.  In  the  sequel,  Fourier  analysis  wiU  be  usrf  o 
analyze  he  cycle  of  an  LDS  of  degree  one  for  constant  coefficient  equations  in  one  dimersMa) 
space  with  periodic  boundary  conditions  when  H  =  2h. 

Let  Eh,  Eh,  be  the  time  marching  operators  on  the  fine  and  coarse  grids,  respectively  Denote 
the  intergrid  transfers  by  jf ,  4.  The  correction  term  r  is  initialized  by, 


=  {if^  Ell^  -  Eh)  vPh 

The  solution  of  the  LDS  of  degree  one  algorithm,  for  n  >  1,  is  given  by, 

'^LDS  =  EhUh  +  {n-  1)E^~^t^ 

=  {Eh  +  (n  -  l)(/f  -  Eh)) 


(6.18) 

(6.19) 

(6.20) 
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Tet  Et.  Eu  i¥J%  denote  the  symbols  of  the  respective  operators.  Let  i/)  denote  the  am- 
pUfication  fkctor  of  the  Fourier  component  6  in  a  cycle  consisting  of  consecutive  steps  on  the 
coarse  grid  before  revisiting  the  fine  grid.  It  is  given  by 

G(0, 1/)  =  II  {Enm  +  [llf{0)El{d)i^{e)  + 

The  amplification  factor  of  any  Fourier  component  should  approximate  as  closely  as  possible  the 

analytic  growth  rate  of  that  component.  In  particular,  when  the  §"7 

in  time,  one  would  like  to  guarantee  that  no  Fourier  component  is  amphfied  by  the  LUb  cycle, 
i.e.,  that  the  amplification  factor  of  the  cycle  A(2/)  satisfies, 

A(u)  =  max  G{9,v')<l  (6.21) 

0e[-7r,7r] 

Moreover,  it  would  be  desirable  if  the  LDS  cycle  was  time-stable  for  any  underlying  discretization 

with  this  property.  .r  i  ?  /ou  m  _  i  a 

For  scalar  equations,  it  can  be  easily  seen  that,  if  \Eh{2vo)\  -  1,  ana 

f\2eo)  =  [iffiOo)  EliOo)  i^{0o)  +  iHiO'o)  iM)  -  Eni^eo)]  M^Oo)  0 


then 


lim  G{6o,  i/)  =  00 


(6.22) 


This  is  in  accordance  with  the  polynomial  error  bounds  derived  in  Section  3.  sho^^i  be  netted 
that  although  the  bounds  were  polynomial,  whenever  condition  (6.21)  does  not  hold,  the  algorithm 
exhibits  a  growth  exponential  in  the  number  of  visits  to  the  fine  grid.  Thus,  the  time-stabihty 
of  the  underlying  discretization  does  not  necessarily  imply  time-stability  of  the  LDS  scheme. 
Moreover,  if  this  condition  is  violated  for  aU  i/  <  i^o,  (e.g.,  when  the  discretization  is  not  dissipative 
enough),  then  for  revisiting  index  u  <  i/o  the  more  frequently  the  fine  grid  is  visited  the  faster  the 

solution  blows  up.  .  x*  x.  -n 

Dissipation  induces  an  exponential  decay  of  the  solution.  Therefore,  a  dissipative  scheme  wiU 

eventually  suppress  any  polynomial  growth  of  the  solution,  i.e.,  for  such  discretizations  A{u)  < 
for  large  enough  i/.  However,  since  the  LDS  algorithm  should  maintain  the  fine  grid  accuracy, 
that  grid  has  to  be  visited  sufficiently  often.  Therefore,  the  discretization  should  have  enoug 
dissipation  to  guarantee  that  condition  (6.21)  holds  for  a  schedule  i/  prescribed  by  the  accuracy 
requirement.  Quite  often,  artificial  dissipation  should  be  added  to  the  coarse  grid  discretization  to 
ensure  the  cycle  is  time-stable.  This  dissipation  may  be  added  either  during  the  time  marching  on 
the  coarse  grid,  or  during  the  initiaUzation  stage  as  weU,  in  which  case  it  should  also  be  introduced 

to  the  fine  grid  scheme.  .  w  j 

The  additional  dissipation  typically  results  in  a  dissipative  cycle.  This  might  hmit  the  method 

applicability  for  truly  long  integration  times.  Hence,  one  should  add  just  enough  dissipation  to 
ensure  that  no  Fourier  component  is  amplified  by  the  cycle.  This,  of  course,  does  not  imply  that  no 
wavenumber  may  grow  during  the  cycle;  to  the  contrary,  imposing  such  a  restrictive  requirement 

leads  to  a  dissipative  cycle  of  limited  practical  interest. 

In  general,  the  length  of  the  cycle  i/  and  the  amount  of  artificial  dissipation  should  be  deter¬ 
mined  simultaneously  to  achieve  the  fine  grid  accuracy. 
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6.3  Staggered  vs  Non-staggered  grids 

Time  dependent  problems  are  often  discretized  and  solved  on  staggered  grids.  The  dispersion 
error  of  discretizations  on  such  grids  is,  typically,  significantly  smaller  than  the  error  in  schemes 
of  the  same  order  on  a  non-staggered  grid.  Hence,  one  can  accurately  solve  on  staggered  grids 
for  substantially  longer  integration  times.  In  the  sequel  the  benefits  and  disadvantages  of  using 
such  grids  when  employing  the  LDS  method  will  be  briefly  discussed.  The  exposition  is  rather 
general  and  is  demonstrated  through  a  particular  example  of  the  linearized  Euler  equation,  which 
has  been  investigated  in  the  present  research. 

A  non-staggered  Cartesian  grid  consists  of  a  set  of  discrete  variables  defined  at  its  vertices  and 
the  system  of  equations  is  discretized  at  these  points.  In  contrast,  on  a  staggered  grid  the  various 
variables  are  located  at  different  positions  in  the  computational  cell  and  the  various  equations 
are  evaluated  at  distinct  points  of  the  cell.  Figure  5  depicts  the  variables  distribution  and  the 
corresponding  computational  cells  for  the  linearized  Euler  equation, 

ut  =  a-Vu  +  cpx 

Vt  =  a-Vv  +  cpy  (6.23) 

Pt  -  a-'Vp+  c(ux  +  Vy) 

where  a  =  (ai,a2)  is  a  two  dimensional  vector.  The  superiority  of  staggered  grid  discretizations 


■  cell  for  V 


cell  for  U 


Cl]  cell  for  P 


Figure  5:  Staggered  grid  discretization  of  the  linearized  Euler  equation 

to  those  on  non-staggered  grid  stems  from  the  fact  that  the  symbol  of  a  general  order  mid¬ 
cell  discretization  of  the  first  derivative  has  the  form  J2T=o  2  ctfe  (where  aj,  depends 

on  the  scheme  coefficients),  while  the  symbol  of  a  central  discretization  of  this  derivative  on 
a  regular  grid  is  {0k  are  scheme  dependent).  The  later  symbol  vanishes  at  tt; 

hence,  poorly  approximates  the  continuous  symbol  iir.  The  phase  error  of  a  non-staggered  central 
discretization  is  larger  than  that  of  a  mid-cell  discretization  of  the  same  order  over  a  major 
fraction  of  the  spectrum  of  wavelength  representable  on  the  grid.  Moreover,  in  non-dissipative 
fuU  discretizations  based  on  central  spatial  discretizations  the  Fourier  components  in  the  upper 
half  of  the  spectrum  have  negative  group  velocity  and  move  in  a  direction  opposite  to  the  physical 
waves.  These  components  can  be  spuriously  excited  by  discrete  boundary  conditions  or  mesh 
refinements  [15,  17].  Therefore,  they  might  deteriorate  the  computation  accuracy  even  when 
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Figure  6:  Refinement  of  a  staggered  grid  discretization  of  the  linearized  Euler  equation  with  zero 
Dirichlet  boundary  condition  in  the  y  direction. 

they  do  not  occur  in  the  physical  problem.  In  contrast,  in  mid-cell  discretizations  all  Fourier 
components  have  positive  group  velocity.  A  considerable  drawback  of  the  staggered  grid  stems 
from  the  variables  placement  at  different  positions  in  the  computational  cell  which  does  not  allow 
implementation  of  characteristic  boundary  conditions. 

The  small  phase  error  of  staggered  grid  discretizations  and  consequently  the  significantly  less 
artificial  dissipation  required  to  suppress  the  polynomial  growth  of  the  IDS  error  give  rise  to 
highly  efficient  algorithms  applicable  for  very  long  integration  time.  The  efficiency  obtained  by 
the  LDS  algorithm  for  the  linear  acoustics  equation  solved  on  staggered  grid  was  25  in  2D  and  66 
in  3D,  (see  Figures  15  and  16). 

The  staggered  grid  has  two  drawbacks  associated  with  the  LDS  implementation.  Both  of  them 
occur  in  the  model  problem  investigated  here  and  are  associated  with  the  intergrid  transfers.  The 
first,  is  that  the  intergrid  transfers  satisfy  the  condition  of  Lemma  2.  Therefore,  for  some  high 
frequencies  there  is  a  0(1)  error  and  enough  dissipation  should  exist  to  eliminate  this  error.  The 
second  problem  is  that  the  solution  interpolation  to  the  fine  grid  requires  extrapolation  in  the 
cells  nearest  to  the  boundary  (see  Figure  6).  This  extrapolation  strongly  amplifies  the  oscillatory 
components  in  the  solution  and  a  way  had  to  be  found  to  circumvent  this  effect  for  our  model 
problem.  Despite  these  deficiencies,  for  the  problems  investigated  in  this  work  the  staggered  grid 
supports  substantially  better  results  than  the  regular  grid. 

7  Numerical  Results 

The  LDS  method  was  introduced  in  the  previous  sections  and  its  properties  were  analyzed.  The 
numerical  examples  presented  in  this  section  aim  at  demonstrating  the  method  strength,  as  well 
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as  supporting  the  claims  and  analysis  performed  in  previous  sections. 

All  the  examples  in  this  work  have  spatial  dimension  two.  However,  generalization  to  higher 
dimensional  problems  is  straightforward  and  was  avoided  due  to  the  heavy  computational  cost  of 
such  simulations. 

In  aU  examples  ^  ^  with  H  =  -^  and  h=  in  Example  6  finer  grids  were  used  as  well. 

The  notation  LDS(d  =  a, 'f  =  b)  used  in  this  section  denotes  an  LDS  algorithm  of  degree  a 
which  revisits  the  fine  grid  after  performing  b  time  steps  on  the  coarse  grid  with  the  LDS  system. 

In  all  relevant  plots,  the  error  is  normalized  with  respect  to  the  solution  norm,  e.g,  the  coarse 
grid  error  is  normalized  by  1^,  where  I^U  denotes  a  restriction  of  the  exact  solution  to 

the  coarse  grid.  Similarly,  the  relative  error  which  measures  the  error  in  an  approximation  to  the 
fine  grid  solution  is  normalized  by  the  latter  solution  norm;  e.g.,  the  coarse  grid  relative  error  is 
normalized  by  where  stands  for  a  restriction  of  the  fine  grid  solution  to  the  coarse 


grid. 


7.1  The  Advection  Equation 

The  first  set  of  examples  consists  of  solving  the  advection  equation  on  the  domain  [0, 1]  x  [0, 1] 
with  periodic  boundary  conditions.  The  equation  is  given  by, 

Ut- a{x,y)ux-b{x,y)uy  =  Q  (7.1) 

In  case  a(r,  t/),  6(x,  y)  are  constant,  an  explicit  solution  of  this  equation  with  initial  data 

u{x,y,{i)  =  U(y{x,y)  (7.2) 

is  given  by 

u{x,y,t)  —  uo{x  +  at^y +  bt)  (7.3) 

Two  instances  of  this  equation  are  solved,  the  constant  coefficient  equation  with, 

a{x,y)  =  1  (7.4) 

b{x,y)  =  0.3  (7.5) 

and  the  variable  coefficient  equation  with, 

a{x,y)  =  1  +  0.3  sin  27ra;  (7.6) 

b{x,y)  =  0.3(1  + 0.4  sin  27rt/)  (7.7) 

In  this  set  of  examples,  except  Example  6,  the  spatial  discretization  was  second  order  upwind  and 
integration  was  performed  by  third  order  Runge-Kutta.  In  aU  these  examples  |  ^  =  0.3. 

In  this  set  of  examples  the  solution  was  restricted  to  the  coarse  grids  by  injection,  and  unless 
otherwise  specified,  quintic  interpolation  was  employed. 

Example  1.  The  LDS  yields  the  fine  grid  solution.  Throughout  this  work,  it  was  claimed 
that  the  LDS  method  yields  the  fine  grid  solution  on  the  coarse  grid.  Clearly,  these  solutions  may 
not  be  identical;  rather,  one  would  like  to  ensure  that  the  error  norm  in  the  LDS  solution  relative 
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to  the  fine  grid  solution  is  similar  to  the  error  norm  in  the  fine  grid  solution  relative  to  the  exact 
solution,  that  is, 

^  \\I^U  -  «'‘ll  ,,  „ 

ll/Ha*!!  ~  11^‘f'll  '  ' 

where  I^U  denotes  a  restriction  of  the  analytic  solution  to  the  fine  grid.  Figure  7  demonstrates 
that  the  LDS  yields  the  fine  grid  solution  in  the  sense  defined  by  (7.8). 

Example  2.  The  LDS  effective  integration  time.  Figure  8  shows  two  examples  of  the  LDS 
applied  to  the  variable  coefficient  advection  equation,  for  smooth  and  oscillatory  data.  It  can  be 
readily  seen  that  the  effective  integration  time  the  discretization  can  be  used  (i.e.,  the  integration 
time  when  the  error  is  small)  drastically  decreases  for  oscillatory  data.  This  figure  should  be  born 
in  mind  as  a  reference,  since  in  some  of  the  next  examples  the  same  equation  with  these  initial 
data  are  used  for  integration  times  significantly  longer  than  the  LDS  effectiveness  time,  in  order 
to  emphasize  difference  between  solutions  for  various  parameters. 

Example  3.  Direct  and  simultaneous  initialization  are  practically  indistinguishable. 

In  Section  .5.1.1  two  initialization  procedures  for  the  correction  terms  were  presented.  These 
methods  vary  in  their  computational  cost  and  memory  requirements.  Figure  9  compares  the 
direct  and  simultaneous  initialization  procedures.  It  can  be  seen  that  at  a  time  longer  than 
the  algorithm  effectiveness  time,  the  LDS  solutions  are  hardly  distinguishable  and  the  observable 
difference  between  these  solutions  is  significantly  smaller  than  the  error  in  the  LDS  approximation 
relative  to  the  fine  grid  solution.  Thus,  it  may  be  concluded  that  these  procedures  have  essentially 
the  same  accuracy. 

Example  4.  The  necessary  and  desirable  orders  of  the  intergrid  transfers.  In  Sec¬ 
tion  6.1  the  orders  of  the  intergrid  transfers  required  to  ensure  the  algorithm  performances  were 
analyzed.  This  analysis  implies  that  for  a  second  order  scheme  the  interpolation  should  be  at 
least  cubic  for  an  LDS  of  general  degree  d.  The  injection  introduces  no  error  to  the  high  fre¬ 
quencies,  thus  is  of  infinite  order.  Figure  10  displays  an  LDS  of  degree  one  and  two  (on  the  left 
and  right,  respectively)  with  linear,  cubic  and  quintic  interpolations.  It  can  be  seen,  for  both 
approximations,  that  linear  interpolation  results  in  an  approximation  worse  than  the  coarse  grid 
solution.  Cubic  interpolation  is  sufficient  to  guarantee  the  LDS  performances,  however,  quintic 
interpolation  yields  significantly  better  results.  This  phenomena  might  be  due  to  the  interpolation 
that  transfers  some  of  the  smooth  components  energy  to  the  fine  grid  high  frequencies  which  are 
strongly  damped  on  the  fine  grid.  This  interpolation  error  is  larger  for  the  cubic  interpolation. 

Example  5.  The  efficiency  of  the  LDS  of  various  degrees.  An  LDS  approximation  of 
higher  degree  provides  a  better  approximation  to  the  fine  grid  solution  than  a  low  degree  one. 
Hence,  it  may  be  used  for  longer  integration  time  before  the  fine  grid  should  be  revisited.  On  the 
other  hand,  for  such  a  scheme  the  initialization  procedure  as  well  as  the  coarse  grid  time-stepping 
are  more  costly. 

Figure  11  addresses  the  question  which  approximation  is  more  cost  effective,  a  low  or  a  high  de¬ 
gree  one.  In  this  figure  LDS  algorithms  of  various  degrees  with  the  same  efficiencies  are  compared 
for  both  smooth  and  nonsmooth  data  and  for  diverse  costs.  The  efficiency  of  LDS(d  —  8), 
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LDS(<i  =  2,7  =  20)  and  LDS(d  =  3,7  =  41)  is  9.4;  while  that  of  LDS(d  =  1,7  =  20)  and 
LDS(d  =  2,7  =  79)  is  15.8.  In  order  for  an  LDS  of  degree  three  to  achieve  efficiency  of  15  the 
fine  grid  should  be  visited  once  in  470  coarse  grid  steps,  while  an  efficiency  of  15.8  can  not  be 
achieved  even  if  the  fine  grid  is  visited  once  in  1000  time  steps.  This  suggests  that  high  efficiency 
can  not  be  achieved  with  high  degree  LDS.  The  efficiency  of  LDS(</  =  3,7  =  350)  is  only  14.6. 
For  short  integration  times  this  efiiciency  can  not  be  achieved  with  LDS  of  degree  three  since  not 
enough  time  steps  are  marched  during  the  simulation,  e.g.,  see  Figure  11  left. 

Figure  11  plots  the  error  in  various  solutions  relative  to  the  fine  grid  solution.  Inspection  of 
Figure  11  reveals  that  for  small  relative  error,  e.g.  ~  0.02  (i.e.,  T  <  1  for  the  smooth  data  example 
and  T  <  0.3  for  the  oscillatory  data),  then  for  smooth  data  the  LDS  of  degree  one  provides  a  better 
approximation  for  both  the  low  and  high  efficiency  scheduling  used.  For  more  oscillatory  data, 
the  degree  two  LDS  yields  better  results  for  the  same  efficiency.  This  might  be  explained  by  the 
fact  that  LDS  of  degree  one  provides  a  fairly  good  approximation  for  the  smoother  components, 
but  not  for  the  more  oscillatory  ones.  The  LDS  of  degree  three  does  not  seem  an  appropriate 
alternative  to  the  lower  degree  approximations. 

It  should  be  noted  that  in  3D  problems,  the  performance  of  LDS  of  degree  two  significantly 
improves.  Thus,  the  efficiency  of  the  abovementioned  scheduling  is  :  for  LDS(d  =  1,7  =  8)  it 
is  13.63  ,  for  LDS(d  =  2,7  =  20)  it  is  15.34  and  for  LDS(d  =  3,7  =  41)  it  is  18.  Moreover,  for 
LDS(d  =  1,7  =  20)  it  is  27.9,  and  for  LDS(c?  =  2,7  =  79)  it  is  38.11.  Thus,  although  in  2D 
problems  there  seems  to  be  little  advantage  to  use  second  degree  LDS  rather  than  first;  in  3D, 
this  changes  drastically.  If  these  schedules  yield  similar  results  for  a  similar  equation  in  3D,  then 
the  LDS  of  degree  two  is  more  cost  effective  for  these  problems  than  the  first  degree  algorithm. 

It  should  be  born  in  mind,  that  the  performances  of  LDS  of  various  degrees  might  change  de¬ 
pending  on  the  equations  or  on  the  data.  Therefore,  it  is  hard  to  give  conclusive  recommendations 
which  method  to  prefer;  and  each  problem  should  be  investigated  separately. 

Example  6.  Richardson  Extrapolation.  In  Section  5.2  the  LDS  method  with  Richardson 
extrapolation  was  introduced.  Although  this  technique  is  limited  in  scope  and,  hence,  should 
be  used  with  great  care;  it  might  be  highly  beneficial  when  it  is  applicable.  Figure  12,  provides 
an  example  when  this  idea  works.  It  consists  of  solutions  of  the  constant  coefficient  advection 
equation  discretized  with  first  order  upwind  forward  Euler  method,  which  is  first  order  in  time 
and  space.  The  equation  is  solved  on  fine  grids  of  128  and  256  points  and  corresponding  coarse 
grids  of  32  and  64  points  on  which  the  LDS  is  solved.  In  this  example  the  normalized  errors  are 
computed  with  respect  to  the  exact  solution.  For  smooth  data,  mesh  refinement  of  the  finest  grid 
by  a  factor  of  two  yields  a  decrease  in  error  by  a  factor  of  1.875  for  the  first  order  scheme  versus  a 
factor  of  5.81  for  the  LDS  with  Richardson  extrapolation.  For  nonsmooth  data  the  fine  grid  error 
is  reduced  by  a  factor  of  1.71  while  for  the  extrapolated  LDS  the  factor  is  4.03.  Thus,  the  error 
reduction  for  the  extrapolated  LDS  is  second  order  both  for  smooth  and  oscillatory  data. 

7.2  The  Linearized  Euler  Equation 

The  next  set  of  examples  involves  solution  of  the  linearized  Euler  equation,  given  by 

Pt  =  a{x,y)p:c  +  K^,y)Py  +  cix,y){u:^+Vy) 

ut  =  a{x,y)ux  +  b{x,y)uy  +  c{x,y)px 
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(7.9) 


Vt  =  a{x,y)V:^  +  b{x,y)Vy-\-c{x,y)py 

This  system  was  solved  for  various  settings.  In  aU  of  them  interpolation  and  restriction  were  of 
sixth  order. 

Example  7.  Staggered  vs.  non-staggered  grid.  Figure  13  demonstrates  the  superiority  of 
the  staggered  grid  discretization  over  the  regular  one  for  the  single  grid  solution  and  consequently 
for  the  IDS  approximation.  It  shows  solutions  of  the  variable  coefficient  linearized  Euler  equation 
on  [0,  Ij  X  [0, 1]  domain  with  periodic  boundary  conditions.  The  discretization  is  second  order 
upwind  for  the  advection  terms  and  second  order  central  for  the  terms  involving  c(x,t/);  the 
integration  was  performed  with  third  order  Runge-Kutta  with  ^  =  0.3.  In  the  fine  non-staggered 
grid  solution  one  can  see  an  oscillatory  component  which  dominates  the  solution.  This  component 
is  visible  due  to  a  large  phase  error  of  the  nonsmooth  wavelengths  on  this  grid.  The  error  in 
these  components  is  even  more  visible  in  the  coarse  grid  and  LDS  solutions.  Clearly,  the  LDS 
is  ineffective  for  this  integration  time.  However,  since  the  error  in  the  non-staggered  fine  grid 
solution  is  already  very  large,  this  is  not  a  real  drawback.  The  large  dispersive  error  is  also 
the  cause  for  the  little  resemblance  between  the  staggered  and  non-staggered  solutions  for  the 
same  data  and  integration  time.  In  contrast,  aU  the  staggered  grid  solutions  do  not  have  those 
oscillations,  and  the  LDS  provides  an  excellent  approximation  to  the  fine  grid  solution. 

Example  8.  The  effect  of  diminishing  CFL.  In  Section  6.1,  it  was  pointed  out  that  when 
^  ^  0,  the  order  of  interpolation  should  increase.  Figure  14  demonstrates  this  phenomenon 
for  the  constant  coefficient  linearized  Euler  equation  on  [0, 1]  x  [0, 1]  with  periodic  boundary 
conditions.  Discretization  on  a  staggered  grid  is  second  order  upwind  for  the  advection  term  and 
second  order  central  for  the  terms  involving  the  c  factor.  Integration  was  done  by  low  storage 
third  order  Runge-Kutta.  It  can  be  seen  that  for  very  short  integration  time  an  oscillatory  error 
prevails  when  ^  =  0.01.  Increasing  by  a  factor  of  10  the  CFL  as  well  as  the  integration  time, 
eliminates  these  oscillations,  yielding  an  excellent  approximation  to  the  fine  grid  solution. 

Example  9.  High  efficiency  LDS  on  periodic  staggered  grid.  The  example  in  Figure  15, 
is  a  solution  of  the  constant  coefficients  acoustics  equation, 

Pt  =  Ux+  Vy 

—  Px  (7.10) 

=  Py 

on  the  domain  [0, 1]  x  [0, 1]  with  periodic  boundary  conditions  discretized  on  a  staggered  grid  with 
the  same  discretization  as  in  the  previous  example,  with  =  0.3. 

This  discretization  has  only  little  dissipation  through  the  Runge-Kutta  scheme.  However, 
since  the  mid-ceU  discretization  provides  an  very  good  approximation  to  the  differential  operator, 
the  coarse  grid  operator  well  approximates  the  fine  grid  operator.  Hence,  only  a  little  artificial 
dissipation  should  be  added.  Sixth  order  artificial  dissipation  was  added  to  aU  equations  of  the 
form  chA^,  with  c  =  0.005.  The  small  dispersive  error  of  the  mid-cell  discretization  enabled  both 
very  long  integration  time,  as  well  as  very  high  LDS  efficiency  of  26  in  2D  and  66  in  3D.  Note 
that  at  this  stage  the  LDS  error  relative  to  the  fine  grid  is  forty  times  smaller  than  the  relative 
error  of  the  coarse  grid. 
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Example  10.  High  efficiency  LDS  on  non-periodic  staggered  grid.  Figure  16  plots  the 
solution  of  the  acoustics  equation  on  [0, 1]  x  [0, 1]  with  periodic  boundary  conditions  in  the  x 
direction  for  all  variables  and  zero  Dirichlet  boundary  conditions  in  the  y  direction  for  u, 

V{x,0,t)  =  V{x,l,t)  =  0  (7.11) 

except  for  this  difference,  all  the  other  parameters  are  identical  to  those  in  the  previous  example. 

Recall  that  high  order  one  sided  interpolation  near  the  boundary  strongly  amplifies  the  oscil¬ 
latory  components.  This  problem  had  to  be  circumvented  in  this  equation  for  the  interpolation  in 
the  y  direction.  The  observation  that  led  to  a  resolution  of  this  difficulty  is  that  for  this  equation 
and  these  boundary  conditions  p,  u  are  symmetric  in  the  y  direction,  and  v  is  asymmetric  in  this 
direction.  These  properties  were  exploited  in  the  design  of  the  intergrid  transfers  as  well  as  in  the 
introduction  of  the  high  order  artificial  dissipation.  The  p,  u  variables  were  symmetricly  extended 
around  the  boundary  in  the  y  direction  and  the  interpolation  and  dissipation  were  calculated  for 
the  extended  solution.  An  assymetric  extension  for  the  v  variable  was  similarly  defined  and  used. 

The  assymetry  v  is  follows  from  the  dual  initial  data  argument,  which  asserts  that  the  boundary 
condition  (7.11)  may  be  viewed  as  requiring  that  a  dual  solution  with  the  same  magnitude  but 
opposite  sign  wiU  constantly  hit  the  boundary  from  the  exterior  of  the  domain  (e.g.,  see  [16]). 
The  symmetry  of  p  in  the  y  direction  follows  from 

Pyt{x,Q,t)=  vtt{x,0,t)-Q  (7.12) 

Thus,  if  initially 

p^(a:,0,0)  =  0  (7.13) 

this  symmetry  is  maintained  for  later  time.  A  similar  argument  holds  for  the  boundary  condition 
at  j/  =  1.  For  the  u  variable, 

Uyt{x  ,  0,  t)  —  Uly{x^  0,  i)  —  Pxyi^X,  0,  /)  =  Pyx{x^  0,  t)  =  0?  0  —  ^  (^•^^) 

Again,  if  the  initial  solution  was  symmetric  in  this  direction,  symmetry  is  preserved.  These 
assumptions  hold  for  the  initial  data  taken  in  our  examples. 

It  can  be  seen  that  the  efficiency  and  accuracy  of  the  LDS  algorithm  were  not  affected  by  the 
imposition  of  non-periodic  boundary  conditions. 

Example  11.  The  effect  of  low  order  artificial  dissipation.  The  parameters  taken  in 
the  example  in  Figure  17  are  identical  to  those  in  Figure  15  except  for  the  use  of  lower  (fourth) 
order  of  dissipation  of  the  form  ehA^,  with  e  =  0.02.  This  choice  of  c  yields  the  same  damping 
of  the  oscillatory  components  as  the  sixth  order  dissipation  used  in  Example  9.  The  resulting 
LDS  solution  does  not  provide  a  satisfactory  approximation  to  the  fine  grid  solution  even  for 
integration  time  significantly  shorter  than  the  one  used  in  Example  9  and  more  frequent  visits 
to  the  fine  grid.  It  can  be  concluded  that  higher  order  dissipation  is  indeed  essential  for  the 
algorithm  performances,  as  lower  order  dissipation  damps  too  strongly  the  smooth  components. 

7.3  The  nonlinear  Euler  equation 

Example  12.  The  LDS  method  for  the  nonlinear  Euler  equation  In  Figure  18,  the 
nonlinear  Euler  equation  is  solved  on  [0, 1]  x  [0, 1]  domain  with  periodic  boundary  conditions. 
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This  system  is  given  by, 


Pt  +  p{Ua;  +  Vy)  +  UP:::  +  V  py  =  0 

Ut  +  {uUj:  +  VUy)+  —pa:  =  0  (7.15) 

C2 

Vt  +  {UV:c  +  VVy)+  —Py  =  0 


where  p  =  p'^ ,  and  7  =  1.4.  Second  order  central  discretization  was  used  for  all  terms, 

with  third  order  low  storage  Runge-Kutta  =  0.3).  Artificial  sixth  order  dissipation  was  added 
with  €  =  0.8.  The  interpolation  and  restriction  were  of  sixth  order  accuracy. 

The  LDS  efficiency  in  this  example  is  17.5  for  2D  problems  and  32.8  for  3D  problems. 


8  Conclusions 

The  Large  Discretization  Step  methods  for  time  dependent  problems  were  presented.  First,  the 
LDS  approximation  was  defined.  It  consists  of  a  system  of  lower  accuracy  discretizations  ap¬ 
proximating  a  more  accurate  time  dependent  discrete  operator.  Error  bounds  on  this  type  of 
approximations  to  linear  problems  were  obtained  for  the  semi-discrete  and  fuUy  discrete  cases. 
These  estimates  hold  for  both  hyperbolic  and  parabolic  equations.  The  research  reported  herein 
aimed  at  deriving  efficient  algorithmic  implementation  of  the  LDS  approximation  for  hyperbolic 
equations,  a  type  of  equations  which  were  not  previously  amenable  to  multigrid  methods.  A 
heuristic  argument  motivated  the  design  of  the  LDS  algorithm  for  nonlinear  problems,  as  well. 

The  LDS  methods  enables  to  obtain  the  fine  grid  accuracy  on  a  coarse  grid  by  adding  correction 
terms  to  the  coarse  grid  equations,  initializing  them  using  the  fine  grid  and  solving  a  system  of 
equations  for  these  terms.  The  accuracy  of  the  correction  terms  deteriorates  at  a  rate  determined 
by  the  coarse  grid  discretization.  However,  since  their  norm  is  significantly  smaller  than  the 
solution  norm,  they  may  be  effectively  used  for  many  coarse  grid  time  steps.  Thereafter,  the  fine 
grid  should  be  revisited  to  compute  new  initial  data  for  them.  Fourier  analysis  was  employed  to 
analyzed  different  aspects  of  the  algorithm;  in  particular,  to  obtain  conditions  on  the  necessary 
orders  of  the  intergrid  transfers. 

The  resulting  algorithm  has  a  typical  efficiency  of  16  for  2D  problems  and  28  for  3D  equations. 
This  efficiency  was  achieved  for  linear  problems  with  periodic  and  Dirichlet  boundary  conditions 
and  the  for  the  nonlinear  Euler  equation  with  periodic  boundary  conditions.  A  particularly  good 
discretization  of  a  linear  equation  yielded  efficiency  of  25  in  2D  and  66  in  3D  problem. 

The  results  presented  in  this  work  are  very  promising  as  for  the  potential  of  the  proposed 
approach  to  tackle  even  more  complex  problems.  StUl,  a  lot  of  research  should  be  carried  out 
to  better  understand  the  method  abilities  and  limitations.  Several  aspects  of  the  LDS  algorithm 
should  be  further  investigated.  The  frequency  the  fine  grid  is  visited  is  important  both  to  the 
algorithm  efficiency  and  to  ensure  the  accuracy  of  the  resulting  solution.  It  would  be  highly  ben¬ 
eficial  if  there  was  a  systematic,  preferably  adaptive,  way  to  determine  when  the  fine  grid  should 
be  visited.  The  dissipativity  of  the  LDS  cycle  is  essential  to  damp  the  polynomially  growing  error. 
For  dissipative  discretizations,  this  requirement  typically  does  not  pose  any  problems.  Quite  of¬ 
ten,  though,  some  artificial  dissipation  should  be  added  to  guarantee  the  algorithm  performances. 
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A  general  way  is  required  to  determine  the  right  amount  of  artificial  dissipation  which  will  sup¬ 
press  the  error  growth  without  affecting  the  algorithm  accuracy.  Boundary  conditions  treatment 
within  the  LDS  does  not  seem  to  pose  particular  problems.  Nevertheless,- a  large  variety  of  such 
conditions  should  be  implemented  and  tested.  In  particular,  a  general  approach  should  be  found 
to  reduce  the  high  interpolation  error  near  the  boundaries. 

The  novel  character  of  these  methods  opens  many  research  directions.  The  application  of  the 
LDS  to  systems  of  conservations  laws  and  shocks  calculations  should  be  investigated.  Another 
interesting  direction  is  to  investigate  possible  generalizations  of  this  method.  In  the  present 
research,  one  solves  for  the  correction  terms  of  a  linear  problem  the  same  equation  as  the  solution. 
This  approximation  does  not  always  yield  the  desired  results,  e.g.,  consider  the  discretization  of 
the  wave  equation  discussed  in  Section  3.  In  a  more  general  setting,  there  might  be  several 
correction  terms  each  satisfying  a  different  equation.  Such  a  generalization  might  significantly 
reduce  the  simplicity  of  the  present  approach;  however,  it  could  be  applicable  to  a  broader  class 
of  equations. 

It  is  expected  that  the  incorporation  of  the  LDS  ideas  into  parabolic  solvers  would  significantly 
improve  their  performances,  as  well.  This  is  another  promising  research  direction. 
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Figure  7:  The  advection  equation,  ut  +  +  0.3  =  0,  with  initial  data  uq  = 

Left  :  (a)  fine  grid  solution  (b)  the  LDS(d  =  1,7  =  20)  solution  (c)  coarse  grid  solution 
Right  :  Relative  errors,  (a)  ||/^u  —  u^||/||/*u||,  (b)  \\u^  —  'u-ids\\I\\'^^\[i 


Figure  8:  Solutions  of  Uj  +  (1  +  0.3  sin(27ra;))Uj;  +  0.3  (1  +  0.4  cos(27rj/))uj,  =  0.  Left  :  Initial  data 
is  uq  =  Right  :  Initial  data  is  uq  =  .  In  both  figures  :  (a)  the  fine  grid 

solution  (b)  the  LDSfd  =  1,7  =  20)  solution  (c)  the  coarse  grid  solution. 


T  =  3.2 


T  =  3.2 


Figure  9;  The  advection  equation,  Uj  +  (1  +  0.3  sin(27ra:))uj,  +  0.3  (1  +  0.4  cos{2-Ky))uy  =  0,  with 
«o  =  Left  :  (a)  the  fine  grid  solution  (b)  the  LDS(d  =  2,7  =  20)  solution  when 

employing  direct  initialization,  (c)  the  coarse  grid  solution  (d)  the  LDS(d  =  2,7  =  20)  solution 
when  employing  simultaneous  initialization.  Right  :  Relative  errors,  (b)  relative  error  of  the 
LDS(d  =  2,7  =  20)  solution  when  employing  direct  initialization,  (c)  relative  error  of  the  coarse 
grid  solution  (d)  relative  error  of  the  LDS(d  =  2,7  =  20)  solution  when  employing  simultaneous 
initialization. 
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T=3.2  T  =  a2 


Figure  10;  Solution  of  Uj+(l+0.3  sin(27ra;))ui;+0.3  (1+0.4  cos(27rj/))uj,  =  0,  with  uq  = 

Left  :  (a)  the  fine  grid  solution  (b)  the  LDS(d  =  1,7  =  20)  solution  when  using  quintic  interpola¬ 
tion  (c)  coarse  grid  solution  (d)  the  LDS(d  =  1,7  =  20)  solution  when  using  cubic  interpolation 
(e)  the  LDS(d  =1,7  =  20)  solution  when  using  linear  interpolation.  Right  ;  (a)  the  fine  grid  solu¬ 
tion  (b)  the  LDS(d  =  2,7  =  20)  solution  when  using  quintic  interpolation  (c)  coarse  grid  solution 
(d)  the  LDS(d  =  2,7  =  20)  solution  when  using  cubic  interpolation  (e)  the  LDS(d  =  2,7  =  20) 
solution  when  using  linear  interpolation. 
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T  =  3.2  T=1.1 


Figure  11:  Relative  error  in  the  solution  of  Uf +  (1  +  0.3  sin(27ra;))u^  + 0.3  (1  +  0.4  cos{2iry))uy  =  0, 
with  Left:  Initial  data  is  uo  =  Right  :  Initial  data  is  uq  =  e-50(2:^+a/^) 

figures  :  (a)  the  LDS(<i  =  1,7  =  8)  solution  (b)  the  LDS(d  =  2,  j  =  20)  solution  (c)  the 
LDS(d  =  3,7  =  41)  solution  (d)  the  LDS(d  =  1,7  =  20)  solution  (e)  the  LDS(c?  =  2,7  =  79) 
solution  (f)  the  LDS(d  =  3,7  =  350)  solution. 
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Figure  12:  Errors  in  solutions  of  ut  +  Ux  +  O.Zuy  =  0.  Left  :  Initial  data  is  «o  =  c“25(a;2+y2)^ 
Right  :  Initial  data  is  «o  =  In  both  figures  :  (a)  error  in  solution  on  a  128  x  128 

points  grid,  (b)  error  in  LDS(d  =  1,7  =  20)  solution  on  a  32  x  32  points  grid  with  Richardson 
extrapolation,  (c)  error  in  solution  on  a  256  x  256  points  grid,  (d)  error  in  LDS(d  =  1,7  =  20) 
solution  on  a  64  X  64  points  grid  with  Richardson  extrapolation. 


Figure  15:  The  acoustics  equation  on  periodic  domain,  with  initial  data  po  = 

Vo  =  0.  In  the  plots  of  the  solutions  :  (a)  the  fine  grid  solution  (b)  the  LDS(c?  =  1,7  =  80)  solution 
when  adding  sixth  order  dissipation  (c)  the  coarse  grid  solution.  In  the  plot  of  the  relative  error  : 
(b)  the  relative  error  in  the  LDS(d  =  1,7  =  80)  solution  (c)  the  relative  error  in  the  coarse  grid 
solution. 
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